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Abstract 

The theory of correlated electron systems is formulated in a form which al- 
lows to use as a reference point an ab initio band structure theory (AIBST). 
The theory is constructed in two steps. As a first step the total Hamiltonian 
is transformed into a correlated form. In order to elucidate the microscopical 
origin of the parameters of the periodical Hubbard- Anderson model (PHAM) 
the terms of the full Hamiltonian which have the operator structure of PHAM 
are separated. It is found that the matrix element of mixing interaction 
includes ion-configuration and number-of-particles dependent contributions 
from the Coulomb interaction. In a second step the diagram technique (DT) 
is developed by means of generalization of the Baym-Kadanoff method for cor- 
related systems. The advantages of the method are that i) a non-orthogonal 
basis can be used, particularly, the one which is generated by AIBST; ii) 
the equations for Green's functions (GFs) for the Fermi- and Bose-type of 
quasi-particles can be formulated in the form of closed system of functional 
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equations. The latter allows to resolve some difficulties of previous versions of 
the DT. Although the expressions for all interactions depend on the overlap 
matrix, it is shown that the theory is formally equivalent to one with orthog- 
onal states with redefined interactions. When the PHAM is treated from the 
atomic-limit side the vertexes are generated by kinematic interactions. The 
latter arise due to non-trivial commutation relations between A-operators 
and come from the mixing, hopping and the overlap of states. The equations 
for GFs are derived in different approximations for vertexes, corresponding to 
the Hubbard-I, mean field, and random phase approximations, respectively. 
The technique is also extended to the case of intersite Coulomb interaction. 
The self-consistent equation for "Hubbard Us" is derived. Both interactions, 
kinematic mixing and hopping, and Coulomb interaction of /- and non-/- 
electrons contribute to the screening of bare U. The diagram techniques for 
pure Hubbard, Anderson models, s — /-model and for iso- and anisotropic 
spin Hamiltonians can be obtained from the developed approach as special 
cases. 

71.25.-d, 71.25,+x, 71.27.+a, 71.28.+d 
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I. INTRODUCTION 



The calculations of properties of real materials in the language of many-body wave func- 
tions are incomparably more effort- and time- consuming than the ones based on the density 
functional theory (DFT). However, DFT is designed for calculation of ground-state prop- 
erties and, by construction, is expected to have limited use if one is interested in dynamic 
properties. The only available regular way to calculate quasiparticle properties is to make 
use of some diagram technique. The latter usually involves expressions for secondary quan- 
tization. Therefore, anyhow one needs some single-electron equation which would generate 
a set of wave functions with proper boundary conditions. The DFT-based ab initio band 
structure calculations (AIBST) seems to be a reasonable starting point at least for the con- 
duction electrons. Therefore, the corresponding Schrodinger equation can be used as the 
one which generates a convenient basis set of wave functions for formulation of secondary 
quantized Hamiltonian. On the other hand it is far from obvious how DFT can be used for 
correlated electrons. Materials which usually are classified as systems with strong electron 
correlations (SEC) show wide spread properties such as magnets with localized and partly 
localized moments, Mott insulators, Kondo lattices, heavy fermion systems, and high-T c 
superconductors. These materials contain elements with an open d— or /—shell which of- 
ten are main source of difficulties in a theoretical description both within the local spin 
density approximation^ (LSDA) and model approaches. The classical example of a ground 
state property which the LDA approach is not able to describe is the Mott insulating state 
where materials with odd number of electrons per primitive cell are insulatorsll. The special 
features which arise in these systems due to SEC are configuration-dependent mixing and 



simpler properties as the atomic-like Hund's rules are not reproduced within DFT. The 
common way to correct the situation in AIBST is to use phenomenological methods like 




spectral weight transferu. Even much 
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Fock approximation: the SIC corrected orbitals (selected at the onset of the calculation) are 
pushed down by the SIC potential, usually to 8-11 eV below Fermi energy, ep, while all other 
/(d)-orbitals remain above ep and are empty. From the point of view of phenomeno logi- 
cal theories of SCES, like the periodical Anderson model (PAM), the picture generated by 
SIC is not satisfactory since usually mixing interaction disappears when /(cf)-level is below 
the bottom of conduction bandit. This makes it impossible to check the PAM mechanisms 
within AIBST. The LDA+U method gives a much smoother picture, since the values of 
the Hubbard U are taken to be 2-6 eV. The models used are also different and there is 
no clear arguments why in one case should us the Hubbard model (i.e. mixing is equal to 
zero) whereas another case the PAM (i.e. hopping is neglected). This raises the question 
whether these arguments can be derived from strict microscopical approach or not as well as 
whether some minimal model can be derived from the exact full Hamiltonian or not. Below 
we will show that when the system contains a subsystem with (quasi)localized electrons, the 
full Hamiltonian H always contains a part, which corresponds to the periodical Hubbard- 
Anderson model (PHAM). Then, representing H in the form H = Hpham + H' we will 
construct a diagram technique which will make possible both calculation of renormalization 
of Hpham by remaining Coulomb terms from H', and to understand to what extent picture 
obtained within the PHAM is accurate. 

Another problem, which is difficult to tackle within the AIBST, is to calculate the spec- 
trum of quasi-particles. Although there are no grounds to consider the eigenvalues of Kohn- 
Sham equations as energies of quasi-particles, this is quite common practice due to its 
relative simplicity. For some cases a simple formula which connects these energies can be 
written0. Part of the dispersion and life-time of quasi-particles, nevertheless, seem to be 
lost. However, the DFT potential treated as a self-energy0 by construction can neither 
provide the energy dependence of the real part of the self-energy, nor reflect the scattering 
of quasi-particles, since it does not have imaginary part. This is especially important for 
the case of SEC, where Hk\(E) sometimes can give three solutions of the equation for the 
quasiparticle energy E = tk\ + ReT,k\(E) (e.g., satellites in the photoelectron spectra of Ni) 
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while the DFT potential, vdft, due to absence of a dependence on energy can provide only 
one solution to this equation. Thus, in order to obtain information about these properties 
of the system one needs a many-body theory. In the present paper we intend to develop 
an approach for the systems with /(d) -electrons in the strong-coupling regime, therefore, 
we choose perturbation theory from the atomic limit. As follows from the above argument 
it is desirable to start with self-consistent AIBST where the localized electrons are treated 
as core electrons. Here we will treat these localized (or almost localized), electrons with 
the help of Hubbard operatorsE'El X pq = \p)(q\. The latter describe transitions from the 
state p of the ion to the state q. The first step then is a formulation of the Hamiltonian of 
the system in terms of these variables. In the case of rare earth metals (Russel- Saunders 
coupling) and a model of homogeneous gas for conduction electrons the work by Irkhin and 
Irkhin§ can be of great help. Particularly, a way to construct many-electron secondary- 
quantized operators, which describe creation and destruction of the /-ion states with the 
help of the standard technique used in atomic spectroscopy, is thoroughly described therein. 
This technique has been also successfully used for derivation of the form of /-density for use 
in DFT-based AIBStH. The case of j — j-coupling is developed much less since an accu- 
rate account of effects of non-local exchange interaction should be performed within a fully 
relativistic field-theory (the case of strong spin-orbit coupling); the cases when a relativistic 
effects can be taken into account via perturbative Hamiltonian, can be treated within our 
approach too. There is one more important motivation for a microscopical derivation of 
the correlated form of Hamiltonians and an effective Hamiltonians. An extension of AIBST 
to thermodynamics is based usually on mapping of the constrained AIBST to some model; 
the parameters of the model are calculated from fitting the model's energy surface to the 
DFT one for some selected configurations of the systemll. Obviously, the choice of the 
model is not unique, i.e., for sufficiently large number of parameters many different models 
can reproduce the DFT-derived dependence of energy on the parameters involved. Besides, 
being within framework of density functional formalism it is very difficult to avoid classi- 
cal approximations, that can introduce uncontrolled errors. Below we shall show that such 



derivation can be performed for arbitrary localized set of orbitals. This means that the pa- 
rameters of the model, as well as correcting terms, are representation-dependent. Probably, 
the best starting step is a DFT-AIBST. As will be seen, our approach for derivation of the 
correlated form of Hamiltonian in the particular case of Russel-Saunders coupling is not fully 
identical to the one, suggested by IrkhirM but their main results can easily be used for our 
purposes. The formalism will be developed in general notations, so that different cases can 
be considered. Compared to the approach developed in Ref.0 the present paper the main 
emphasis will be put on the next logic step of development of the theory, namely, on the 
derivation of a regular diagram technique from atomic limit for the evaluation of Green's 
functions within the basis chosen. 

The fact, that SEC have local, site- centered nature dictates a requirement of the ap- 
proach: the strong on-site Coulomb interactions should be taken into account first. This 
can be done much easier within the representation of some site-centered wave functions. 
Then the question arises immediately: should one first orthogonalize the wave functions 
centered on different atoms and then take into account correlations or, vice versa, first to 
take into account local interactions and then treat somehow non-orthogonality? The answer, 
of course, depends on the ratio of the bandwidth, generated not only by non-orthogonality, 
but by all possible mechanisms (due to self-consistency) and the largest matrix element of 
local Coulomb interaction. It is clear that a big enough Coulomb interaction can a cause 
strong shift of the levels of orbitals, that, in turn, can strongly change the overlap with 
neighboring atoms. This picture has been many times confirmed by the calculations, which 
use LDA+U-, SIC (self-interaction correction), or Hartree-Fock approximation: the orbital, 
which was delocalized before application of, say, SIC, becomes completely localized after 
applying SIC, loosing simultaneously all mixing and overlap with neighbors. Thus, for such 
materials the answer seems to be clear and conceptually important: first should be taken 
into account SEC and only after that the overlap, mixing and hopping should be treated. 
Strictly speaking, the final band structure should not depend on our choice of perturbation 
theory: from atomic or itinerant limit we develop perturbation theory. In the latter case a 
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necessity arises to generalize the three-body Faddeev equationsE-J to the rz-body ones and to 
work with many-electron GFs. Although this is not done yet, one can expect that an analysis 
of the contributions from the closed channels of scattering should raise the same problems, 
which we have deal with when start from the atomic limit. Returning to our approach, we 
immediately find that the non-orthogonality of orbitals causes a problem: a regular dia- 
grammatic approach for the description of correlations of electrons in open many-electron 
inraatomic shells, for the case when the system is near the atomic limit and the basis used 
is non-orthogonal, has not been developed yet. We have not found a method that allows to 
treat the non-orthogonality within the versions of the DT, developed by Kuramotoi! and 
GreweS as well as for other version of the DT, that has been developed by Westwansky and 
PawlikowskiEl and by Zaitsev§ (WPZ) for the s-band Hubbard model (the group SU(4)). 
The problem is that in the case of a non-orthogonal basis the picture of interaction cannot 
be used since we are not able to calculate the time dependence of the conduction-electron 
operators Ck\(t) and of the X n (t) -operators even for the zero Hamiltonian (18): c(t) and X(i) 
are coupled by the overlap matrix. In the WPZ-technique one more open question exists. 
Namely, as has been shownil, compared to the conventional DT for fermions and bosons, 
the set of diagrams in the DT for the s-band Hubbard model is not unique, being dependent 
on the order in which different X-operators are picked up from the chronologically ordered 
product (T{X 1 (r 1 )...X m (r)..X n (T n )}) for Wick's decoupling. In order to avoid under- or 
double counting, this way of the formulation of the DT must be complemented by a system 
of priority, or hierarchy, for the X-operators. The recipe for the s-band Hubbard model is 
developed. However, it is difficult to use this recipe in practical calculations for the d- and 
/-systems due to very large number of the electron transitions involved (this number quickly 
grows with the number of electrons in the shell). Another important moment is connected 
with existence of closed system of equations for GFs. The latter give us an opportunity to 
introduce the accurate definitions for self-energy and vertexes. We shall show below that an 
absence of such definitions in earlier theories can lead to erroneous results. As will be seen, 
the mentioned difficulties can be resolved by making use of a properly generalized method 



of Baym and KadanoffS0 where the diagram technique can be constructed via equations 
of motion. Within this technique the hierarchy is automatically established by the choice of 
the variables, Grassman or boson ones, which are used in the Hamiltonian of external fields. 

The aims of the present paper can be formulated as follows: 
1) To derive an approximate Hamiltonian in a local many-electron representation, which 
allows to use as a reference point the LDA Hamiltonian, Hlda, for a description of delo- 
calized electrons, but the many-electron description of the localized or, almost localized, 
electrons. To separate from the total Hamiltonian the one corresponding to the generalized 
Hubbard-Anderson model in order to elucidate the microscopical origin of the parameters 
of this model. 2) To derive a closed system of equations for the Green's functions and to 
develop a diagram technique, that allows one to work with a non-orthogonal basis set which 
is free from the problem of a hierarchy between Hubbard operators. 

3) To derive a few " common" approximations for the case of a generalized Hubbard- Anderson 
model within a non-orthogonal basis. Particularly, to derive the self-consistent equations 
for the set of "Hubbard U-s. 

4) To inspect the role of remaining terms of the Coulomb Hamiltonian in formation of the 
Hubbard- Anderson model and to make an attempt to formulate "minimal" model, which 
describes formation of the low-energy physics in the SEC system. 

5) We will finally compare our equations with some of available results of other authors. 

It is worth noting that i) although our derivation of perturbation theory is based on the 
ideas similar to ones used by Ruckenstein and Schmitt- Rinkil (Schwinger's source theory!!) 
the resulting Green's function technique is closer to the one used by ZaitsevS and Izyumov 
et al.H (within "boson" hierarchy) since we are using local Green's functions as a zero-order 
limit instead of the itinerant limit used in ref. 240; ii) we suggest to use the self-consistent 
LDA potential for generating the set of wave functions and eigenvalues as a convenient first 
step, taking in a proper way symmetry and detailed charge distribution of the material of 
interest, but not the conceptional framework of the density functional theory (DFT). We 
emphasize that we will use a more or less standard many-body approach, which is different 
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from DFT ideology: in our case the expressions for the ground state energy, free energy and 
excitations should be found directly from the calculated GFs. 

The content of the next sections is the following. In Sec.IIA, assuming that the self- 
consistent solution of the band problem within LDA or some other convenient potential with 
the /(d) -electrons treated as core is found, we use a local set of site-centered non-orthogonal 
functions for a transformation of the total Hamiltonian to a form that is convenient for 
our further purposes with representation of non-orthogonal field operators. Then we show 
that the Green's function, defined on these non-orthogonal operators, reproduces the LDA 
problem. In section IIB we introduce the many-electron representation and we derive the 
generalized Hubbard-Anderson model corresponding to the secular equation of the LDA 
problem plus those terms of Coulomb interaction which are responsible for formation of the 
/(d) /shell of the ion. 

Sec. Ill is devoted to the diagram technique itself. We start with a derivation of commu- 
tation relations between many-electron operators of the /-shell and delocalized electrons. 
Then with the help of these relations the equations of motion for the operators of quasi- 
fermion and quasi-boson intra-atomic transitions are derived. The latter equations allow 
us to find out what type of Schwinger's sources should be introduced into the theory and 
to derive a closed system of equations for Green's functions (GFs) in terms of functional 
derivatives of generating functional. 

In Sec. IV we rewrite the theory in a more compact form, introducing the self-energies 
and vertexes, and derive a few approximations. The equation for the self-energy is written 
here in terms of functional derivatives, which is convenient for iterations. After that the 
equations for GFs within the " Hubbard- F approximation are formulated. This approxima- 
tion corresponds to the vertex equal to zero. Next approximations require knowledge of the 
vertex. Therefore, we formulate the equation for the self-energy in terms of the vertex. This 
allows us to derive a mean field theory, beyond the Hubbard-I approximation, in which the 
energies of /-transitions are self-consistently renormalized by kinematic interactions caused 
by mixing interaction and hopping. The abovementioned self-consistent equation for Hub- 
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bard U-s is composed of these equations. At last, by making use of transformation to a 
description in terms of an effective field, we formulate the equations for the " random phase 
approximation" (RPA), which describe the self-dressing of kinematic interactions. 

In Sec. V we consider some of remaining terms of Coulomb interaction. The purpose of 
this section is two-fold: i) to generalize the equations and method for the case of Coulomb 
interaction; ii) to discuss the correction from Coulomb interaction of /-electrons and conduc- 
tion electrons to the self-consistent equation for "Hubbard Us" . The paper is completed by 
Sec. VI where discussions and conclusions are given. Applications of the presented technique 
to concrete problems will be given in forthcoming papers^. 



II. HAMILTON! AN AND THE HUBBARD- ANDERSON MODEL 

A. Non-orthogonal representation 

In order to have a possibility to correct the band structure and the ground state proper- 
ties, obtained by the LDA calculation in a regular way we first introduce the exact Hamil- 
tonian of the system in second quantization using some complete set of wave functions, for 
example, plane waves. One can expect that the many-body corrections to the spectrum and 
wave functions of conduction electrons will be smallBH, which motivates to start evalua- 
tions with the LDA Hamiltonian. We note, however, that the construction described below 
allows to start from any convenient single-electron potential (e.g. a self-interaction corrected 
one). Let us separate the LDA part of full Hamiltonian and develop a many-particle per- 
turbation theory for the Green's functions over the deviation between the exact and LDA 
Hamiltonians. 

The corrections which we are going to discuss come mainly from on-site-localized 
Coulomb repulsion in the open d- or /-shell. From a technical point of view a serious 
disadvantage of local representation is the non-orthogonality of wave functions, centered 
on different sites. The corresponding diagram technique should necessarily include over- 
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lap matrices. The same reason, however, turns out to be an advantage since it allows to 
easily separate the on-site interactions from the contributions coming from other sites. In 
the present paper we restrict ourselves with a non-relativistic consideration, i.e., the total 
Hamiltonian of the system we shall start with is: 



H = (T + H en ) + H ee + H nn 

,2m Y /• hj 



+ ^EE/^/ ^(r)^(r')-^-^(r')^(r) + I £ (1) 

o- cr' II I i i'l 



where the first term includes the kinetic energy and the electron-nuclei interaction, the sec- 
ond term the electron-electron interaction and the last term the nucleus-nucleus interaction. 
It should be emphasized that the electron operators ip a ( r ) are constructed from some com- 
plete set of functions and, therefore, the field operator and, correspondently, Hamiltonian 
H, are exact. Let us now introduce the set of wave functions, (fik\( r ), which are the solution 
to the LDA band problem: 

p 2 

^2m + VLDA ^'^ kX ^ = ^x^kxir), (2) 

where vlda(j') is the LDA effective potential, e^A the Kohn-Sham eigenvalue, and A is a 
band index. We will use below the solutions of this equation only as an almost complete 
set, on the one hand, in order to reformulate our full problem into convenient for our 
aims representation , and in order to express the energies of conduction-electron quasi- 
particles in terms of the Kohn-Sham eigenvalues, on the other hand. Since the set <pk\{ r ) is 
orthonormalized, corresponding annihilation operators can be obtained simply by projection 
from the full field operator, ak\ a = j dr<pl x (r)ip a (r). Next, we separate the full field operator 
into two parts, 

Mr) = ^ DA \r) + (Mr) - i>? DA \r)) = <Mr) + ^(r), (3) 

where 
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^{r)=^ DA \r) = Y J <Pkx{r)a kX . (4) 

k\ 

Let us now insert ip = ip\ + ip2 into the exact Hamiltonian (1), and of all terms which arise 
we take into account only the ones which contain only We stress that the approximation 
made here should be kept in mind, since for a description of some experiments, for instance 
the ones involving large-energy perturbation, (e.g. photoelectron spectroscopies) this set 
may be insufficient. For these cases a consideration within this reduced Hamiltonian would 
miss essential physics of the phenomena (essential additional terms in the Hamiltonian). 

The obvious, strict way to construct the standard weak-coupling perturbation theory 
which takes into account the effects of Hubbard U and starts from H LDA as the zero Hamil- 
tonian is to represent H in the form: 

H = (T + H en + V LD a) + H nn + (H ee — V LDA ) = H + H int = H + H int . (5) 

It worth noting that although we use Vlda one can start with any convenient potential. This 
opportunity may be especially important for consideration of problems with constrained 
boundary conditions (like, for example, quantum dots). In consideration below we will 
work, however, with the approximate Hamiltonian, having neglected all the t/> 2 - dependent 
terms: 

H = J dr4*(r) + v LDA (r)^j ^(r) 

+ (H nn + H ne + H ee - J dr^[ a (r-)v LDA {r-)4> la (r-)^ , (6) 

with conventional notations. 

The further steps are the following. The first one is to construct expansion of field 
operators in the formalism of secondary quantization in such a way that the equation of 
motion for the Matsubara (or, Baym-Kadanoff) GF in the XjL- re P resen t a ti° n reproduced 
Eq.10 automatically. The latter means that we have to define the field operators in a proper 
way, then, using this definition, to calculate the commutation relations between them and 
rewrite the full secondary quantized Hamiltonian in this representation. Then we will use 
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these fermion operators for constructing the many-electron operators, which describe SEC, 
while the elementary commutation relations for these fermions will be used for calculation 
of the commutation relations between different complex combinations of them, which arise 
in the theory. In the second step we will separate the Hubbard-Anderson model from the 
full Hamiltonian. This should help to clarify the microscopical origin of the parameters 
of these models, and the nature of remaining Coulomb interactions. In third step we will 
generalize the Baym-Kadanoff formalism to the case of these complex non-Fermi and non- 
Bose operators which are generated by the equations of motion. 

In cellular methods of AIBST the Bloch functions (pk\(r) of the band A are often repre- 
sented in the form of expansion over site-centered functions Xjhij') = Xl{t — Rj)- Here L is 
complex index: L = {I, mi a}. The secular equation in this representation has the form: 

(EO jLJlLl - h jUhLx )x h L x (r) = 0; (7) 

here J x J1 ^ 1 = (XjxIXjiLi) is overlap matrix. The functions if and x are connected by 
relation 

<Pkx(r) = £ u* a \k)Z2(k)e ikR i XL (r - R 3 ) ee 5>fec,- L . (8) 

a jL jL 

The coefficients Z here come from the Choleski's factorization of the overlap matrix, O = 
ZZ, while the Hamiltonian h of corresponding secular problem, 

h aa '(k)(f>ka'(r) = E a (k)(f) ka (r) (9) 

is connected with the matrix element h as follows: 

haa'(k) = Y,lZ~\k)] a L h LL ,(k)iZ-\k)} a L :. (10) 

LV 

At last, the equation for corresponding to the Eq.7 fermion GF F has the form: 

(EO — h)jL,j"L"Fj"L",j'L> = 5jj'5LL>- (11) 

Having known the expression for if it is easy to find definition for a^. Let us insert Eq.8 
into expansion for the field operator (4): 
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iMO = EtE < X (k)Z2(k)e ikR i XL (r - Rj)]a kX = 

kX a jL 

Y,X jL [Y:K\k)Z2(k)e ikR ia kX } E^ Wt , (12) 

jL k\ jL 

i.e., the operator of destruction of electron on the orbital (jL) should be defined as the 
following combination of normal band operators: 



a 3L = Y, U k L \ a k\- (13) 

fcA 



It follows from this definition that the anticommutator of these operators is equal to the 
inverse of the overlap matrix: 

{a jL ,a], L ,} = (O-^jLju, (14) 

OjL,j'L> = (XjL\Xj'L'}- 

Since in practical calculations only finite number of the functions XjL can be used, it is 
clear that the accuracy of this approximation is determined by the possibility to represent 
(^-function in terms of the functions x, namely, 

5(x-x')& E XjL(x)0-l jlL ,x*> L >(x'), 

jl,j'L> 

The terms containing ip2 can alway be taken into account additionally. From a technical 
point of view the case of /-electrons differs from the case of ci-electrons only by the number 
of correlated orbitals involved, at least in the initial steps. Let us now discuss the case of 
/-electrons and use for the operator with L = {/ = 3, mi, a} the notation fj mia , whereas 
all the non-/'s are called c jL , a jL = 5 L ^ mia) f mia + (1 - 5 L ^ mia ))c jL ) . Then the LDA part 
of the Hamiltonian can be written in the form: 

HlMTO = E e k\a<A\a a k\a = E HjL,j' L'c] L Cji L i + E e °ffjmicrfjmi<T 
kX j,L,j',L' j,mi,a 

+ E tjmhj'm^fjmtcrfj'mft 
jmicrj'm[ 

+ E [Hjfyj'<r c ]Lfj'm t a + Hj mia .j, L fj mi(T CjiL\- (15) 
j,L,j'mi,a 
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Here the hopping matrix element and the position of the center of the /-band are given byEil 

k 

4 = Y, H m l a,mA k )- ( 16 ) 

k 

In the case when in ab initio band structure calculations (AIBST) the /-electrons have been 
treated as core electrons, eg is just single-electron level in the "LDA"-atom, since in this 
case hopping and mixing are set equal to zero. Clearly, we are not allowed to diagonalize 
this Hamiltonian directly, since we deal with a non-orthogonal basis set. 

Let us now write down the equation of motion for the Matsubara's Green's function (GF) 

G jL ji/{iw) = -i(r[a iL (t)a] L (t')])^, (17) 

in this representation. We find that in the zero order approximation this equation has the 
desired form 

[lUJ n OjL,j'L' — HjLj'L^GjLj'L'tyuJn) = d~jj>d~L,L', (18) 

since the expression in the square brackets exactly reproduces the matrix of the LDA secular 
problem. 

In order to avoid confusion we remind that, as follows from Eqs. (11?) and (14?), a 
straightforward diagonalization of the matrix Hji^iy would lead to eigenvalues which do 
not coincide with the band-structure ones since the functions XjL are not orthogonal to 
each other, (jL\j'U) = Ojlj'l 1 7^ &jjio~LL> (0> an d the overlap matrix O is hidden in the 
commutation relations (14). 

The other warning is the following: although the poles of the GF in Eq.ll coincide with 
the eigenvalues of the Kohn-Sham secular problem (what we want to have here), one cannot 
obtain the total energy E tot from the standard formula via GF and its self-energy. Therefore, 
the LDA potential in the equation for GF cannot be considered as self-energy. 

The density of charge p a (x) which enters the H LD a, can be found, as usual, from GF 
(10), 
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p a (x) = (iplix^ix)) S3 (lpl a (x)lp 1>a (x)) = lim XjMXj'L'(x)G jL ,j'L'(t,t'). 

jLJ'L' 

As seen from (5-6) the perturbation theory from weak-coupling limit can be constructed 
with respect to deviations (H — Hlda)- In the case of expansion near atomic limit this is 
not possible. 

B. Many-electron representation and periodical Anderson model. 

A guiding idea on how to introduce this representation can be obtained from the known 
considerations of the s-band Hubbard model. Namely, the decomposition of the fermion 
operator cj = c-f • (1 — njj + cj • n^, is produced by the Hubbard repulsion which makes 
energies of the transitions Cf(l — nj.) and c^h^ different. The physical meaning of this 
decomposition is that every electron should "know" about populations of other orbitals at 
the same site. It is clear from this speculation that, as an example, for an /-electron in the 
orbital mi = — 3, cr =| we have to take into account in the corresponding decomposition as 
sum of all possible products: 

/-3,T = /-3,t ■ (! - ™-2,t) • (! - ra_i lT ) ■ ••• • (1 - n +3tl ) 
+ f-3,1 ■ n_ 2 ,T • (1 - n_i,t) • ... • (1 - n +3jl ) 
+ 

+ /-s,T ■ n-2,t ■ ■ ••• ■ "+34 ( 19 ) 

This transformation is, obviously, identical. In the next step, inserting this into the total 
Hamiltonian we take into account correlations by making use of the properties of the Fermi 
operators: f 2 = (f^) 2 = 0, h 2 = n. The sums of the diagonal terms in the Hamiltonian, 
which have the same operator structure, give energies of corresponding many-electron states 
in the zero approximation. This step is still exact. However after this transformation 
returning to the form (6) and (15) of the Hamiltonian becomes somewhat difficult. The 
fermion GF in the atomic limit now acquires the form0: 
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y**)(iu,) = (Tf v fi)\ 



E 



(|r w |/,|r ra+1 )| 2 (iy rn + iv r „ +1 ) 

iu - (E Tn+1 - E Tn ) 



(20) 



where the spectral weights, say, for (i-electrons, contain non-decoupled correlators, like 



and so on. Here E^ n is energy of n-electron configuration r„ of d(/)-ion. The transi- 
tions, described by different terms in this expansion, have different energies, and the energy 
separation between different terms is often so large, that their contribution to almost all ob- 
servables is negligible. These terms can be neglected without changing the physics. The new 
reduced Hamiltonian can still be treated within a single-electron approximation, if we put 
all population numbers of the /-orbitals equal to integer numbers (0 or 1). The procedure, 
described above, actually, gives a formal derivation for the approach, exploited usually in 
band structure calculations for materials with well-localized 4 /-electrons : /-electrons are 
placed in the core and their occupation number is taken to be equal to an integer number. 
For d-s and 5/-s, and even for the first part of 4/-series this procedure is often not sufficiently 
accurate. Thus, if either mixing or hopping are allowed, or some of the many-electron states 
have close energies, we come to the case of strong electron correlations. The key question 
is, of course, how large is Hubbard U in a particular compound with d- or /-electrons under 
interest. It should be noted that the definition of this magnitude and the way of calculation 
of it in different papers is different Hubbard in his first paper0, i.e. just the 

matrix element, which stands in the Hamiltonian. Note, however, that if it is convenient, 
any "preliminary" value can be chosen for H , since it will be taken into account later via 
H — Hq. Sometimes the shift U — > U mo dei is desirable, since this decreases the strength 
of perturbation. Since we started from an exact Hamiltonian, we have an opportunity to 
calculate U mo d e i in terms of Coulomb interactions, existing in the system. We will return to 
this question after including Coulomb interaction into our diagram technique. 

As we noted above, an alternative way to transform the Hamiltonian to the many- 
electron form, which describes the correlated motion of /-electrons in the case of Russel- 
Saunders type of coupling used by Irkhin0. After the transformation is done, the terms 



iV r , 



((1 - n0(l - n 2 )...(l - n v -\)n v {l - n v+1 )...{l - h w )) 



(21) 
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in the expansion, which clearly give small contribution can be neglected (for example, one 
can hardly expect that the Ni ion can be found in 5— or 6— valent states). Below we will 
use an equivalent, but more convenient way, compared to the above described method, to 
introduce the many-electron representation. 

In the following part of this section we shall perform the following three steps: 1. In- 
troduce a many-orbital representation. 2. Rewrite our Hamiltonian in this representation. 
3. Separate out the main terms of the perturbation Hamiltonian. These steps will allow 
us to derive microscopically the multi-orbital periodic Hubbard-Anderson model for this 
particular representation. 

Let us now construct the orbital representation for atomic states. It is convenient to use 
the following short notations for the /-orbitals. The orbitals with spin "up" starting from 
mi = —3 to mi = 3 we shall label by /i — 1,2,..., 7, whereas the ones with spin "down" 
from mi = —3 to mi = 3 by fi = 8, 14. For instance, \m\ = —2, j) = |/i = 9). The states 
of an /-ion with different number of electrons in the /-shell will be described in the orbital 
representation as |r ) = |0), |r M ) = |//) = /t|0), \r^) = = (^/J/t|0) + 

and so on. Here 9 Ufl — 1 if v < // and 9 Ufl — if // < v. In other words, the following 
convention about the order of orbitals in any many-particle state \T n ) = /i2, /i„) = 
/1/U2---/LI0) is taken: /ii < /i 2 < ... < /i n always holds. The advantage of the orbital 
representation is that it presents the simplest possible way to transform the Hamiltonian to 
the atomic representation and gives a very simple form for zero fermionic Green's functions. 
Another advantage is that almost all rare earths have single-determinant ground state. On 
the other hand, the set of wave functions |T n ) is a complete set for any fixed n, therefore, 
more complicated states can be easily combined of these wave functions. This representation 
can be good starting point for description of ci-electrons in 3<i-materials or for 5/-electrons 
in actinides. For the case of 4/-electrons the Russel-Saunders coupling should be used, i.e. 
many-electron states have to be constructed according to the scheme (see recent review by 
Irkhin and Irkhin0 and references therein): 
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r»> = 4 jo) = 4= £ G^C^JlAljO) (22) 



n r 

where Gr™^ = G s "^™^™_ 1Qin _ 1 are fractional parentage coefficients (a is Rakah's seniority), 
and Cp" , = Cf n ' M M w^Cc"' M ?, 1/9 are Clebsch-Gordan coefficients. For actinides LS'- 

1 n-l in-lMn-l™ b n — lll n — 1,1/2 

coupling should be used (with spin-orbit coupling and other relativistic terms added). For 
3<i-materials, probably, crystal-field combinations are a proper choice. Thus, since different 
materials require different combinations of wave functions (or, operators), below we will use 
general notation for expansion of /-operators, so that a combination needed for diagonal- 
ization of Ha. 



ton 1 



H = V e° a* a ■ 

+ 2 ^MlM2M3M4 a j/Xi a ]^ 2 a iA t 3 a iM4 — E U MlM2 a iw a j>2> (^3) 

is assumed to be used. Here the e M are core and /-levels from LDA-calculation. Note, that 
in the case of Russel- Saunders coupling the Coulomb part of interaction can be rewritten 
in terms of Slater integrals and conveniently separated into the terms describing spin-spin, 
orbital-orbital interactions, etc.0. 

Let us now introduce a description in terms of the Hubbard's operators X rr ' = |r)(r'|. 
In order to see explicitly the process involved in the particular term of perturbation the- 
ory, it is convenient to use different notations for the diagonal and non-diagonal opera- 
tors, distinguishing also quasi-fermionic and quasi-bosonic operators: For diagonal opera- 
tors we use the notation: h T = |r)(r|. For non-diagonal Fermi-like we use the notation: 
X r ™' r ™±i = |r n )(r^ ±1 |. For non-diagonal Bose-like, with no change in the number of parti- 
cles, we use Z: Z Tn,v ' n = |r n )(rjj and for Bose-like operators where the number of particles 
changes by an even number we use Z Vn,Vn±2 = \T n )(r n ± 2 \- Below we shall use also shorter 
notations for the transitions; for the Fermi-like transitions we denote X a with a = [T n , r' ra+1 ] 
, and for Bose-like transitions we use capital as super scripts: Z^ with £ = \T n , T' n ] (includ- 
ing h r = Z T ' r ), or A = [r n ,r^ ±2 ] . The transitions in opposite direction will be denoted 
by a bar: a = [T' n+1 ,T n ]. We also introduce special notations for the anticommutator of 
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AT-operators and its average: 

TrAe~P H 

{XIX)} = Qf and ({Xf,X*})=Qf with (A) = — ^ . (24) 

Here /3 = 1/T, and T is the temperature in units ks = 1- The fermion operator can be 
written in this representation as 

/ i7 = u,) a x? = E /o^X' 71 + E / 7 V*7' r + E A + - , (25) 



71 



where 7 describes single-electron states, T and A stand for two- and three-particle states cor- 
respondingly. The matrix elements /^ 7l = (0|/ 7 |7i), /J r = (V|/ 7 |r), / r Y , A = (r'|/ 7 |A) of the 
fermion operator / 7 (using many-electron states \A) = \T), |A)) describe selection rules 
and reflect the representation chosen for description of /(d)-states: j — j-, Russel-Saunders 
coupling or the orbital one. In the case when the /-(non-/) overlap integrals are neglected, 
these coefficients acquire the simplest form in the orbital representation. For example, 
(7 ; |/ 7 |r) 7^ only if |T) = |7'7). Therefore, in the expansion (25) the sum over double- 
electron states T is actually, X 1 ' ,<k1 ' ^{1—5^ n ) . In the last term in (Eq.25) the summation 
over three-particle states A includes only the states with A = (r", 7). Thus, the sum over r" 
contains only those of states T', which do not enter the second term of Eq.(25). For the rare 
earth ions the /-electron coefficients of expansion over X-operators can be expressed via 
Clebsch-Gordan and fractional parentage coefncients@0, (y1) a =[ r n-i> r "] = y^nG^^C^^ 7 - 
However, when either /-orbitals of different sites, or /- and (non-/)-orbitals, are not or- 
thogonal to each other, these coefficients contain inversed overlap matrices. The Hubbard's 
representation, taking into account correlations, builds into the theory the physical fact that 
in a system with SEC around atomic limit, the single electron transition, described by a 
fermion operator, does not exist as a simple one-electron excitation. Instead it is split into 
many intra-ion transitions. Formally this expansion includes the transitions between all 
states with (n + l)and n electrons. However, as was pointed out by Hubbardll, the terms 
with large deviation of the population number rif from its equilibrium value are separated 
by a very large energy gap from the ground state and therefore, if they even exist at all, 

20 



they do not give contribution to the formation of a crystal (and even liquid). 
Let us change every /-operator in the total Hamiltonian by its expression in terms of X- 
operators (Eq.25) and separate the intra atomic part. The terms of Hamiltonian which 
contain more than one Hubbard operator, belonging to the same ion, should be simplified 
according to the multiplication rule for Hubbard's operators, X pq X r - s = u qr Xj S , what takes 
strong interactions into account algebraically (here v qr = (q\r) is sum determinants consist- 
ing of matrix elements of corresponding to these states inversed overlap matrices; v qr ^ S qr 
if 0~ l 7^ I . After that we have to collect the terms which have the same operator structure. 
Performing this transformation with the LDA-part of the Hamiltonian (6), and adding the 
/ — /-Coulomb single-site interaction, we come to the first-step" Hubbard- Anderson model: 

Wo = £^r fc f+ £ H jLJfL J jL c fLf , (26) 

'•I jL,j'L> 

n™ AM = £ t%X]X), + £ W jL , f J jL X« + H.C.. (27) 
The matrix elements of the Hamiltonian are discussed in Appendix A. 

We need also the expression for the total number of electrons in this representation: 

N e = Oj 1 L 1 ,j 2 L 2 c ] 1 L 1 c j2L 2 + OjiL 1 ,j2M 2 (.fM 2 ) ac ] 1 L 1 Xj 2 

+ O jl M 1 ,j 2 L 2 Uii 1 ) aX h c hL2 + (^M)r,r^J 
+ (! - ^hh){fM 1 ) a OhM 1 ,j 2 M 2 {fM 2 ) b X^ i X b j2 . (28) 

As seen in the equations above, due to the non-orthogonality between states, part of the 
electrons are in mixed states. 

Now our Hamiltonian is fully defined, i.e. the operator structure is derived and a recipe 
for a calculation of all matrix elements has been given. This allows us to continue to the 
next step, i.e. to the construction of the diagram technique. The reason why we call the 
Hamiltonian (26-27) as the "first-step" one is that remaining terms of Coulomb interaction 
renormalize these parameters. We will return to this question later. First we will consider 
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the Hubbard-Anderson Hamiltonian. As we said, due to very large number of terms in 
full Hamiltonian we will construct first the diagram technique for the Hubbard-Anderson 
Hamiltonian, Eqs. (26), (27). After that we will return to Coulomb interaction. 

III. DIAGRAM TECHNIQUE 

A. Commutation relations 

Since our aim here is to calculate the spectra and spectral weights of the quasi-particles 
by means of Green's functions, the commutation relations between the operators involved 
is a necessary ingredient. Let us start from the fermion-like excitations and derive the 
equations of motion id t A = [A,H] for the operators A = Cji,XJ, for the case when these 
operators for different sites neither anticommute with each other nor commute with hr- 
The correlations here manifest themselves via commutation relations on the site, to which 
the /-electrons under discussion belong. The results of a commutation [XJ, 1 , h^,] , (which 
gives two Fermi-like X-operators), and an anticommutation {X,X^}, (which gives Bose-like 
operators h, Z) are properties of the particular algebra used. It is convenient to have this 
information, written in terms of the structure constants of the correspondent algebra: 

{X a ,X 1 } = efh r + efZ A + efZ A = £ efz*, (29) 

Z=t,a,a 

[X a ,h r }=e a a fX ai . (30) 

When the interaction between X- and c-subsytems is switched off, in this particular (orbital) 
representation the constants e can have values 0, 1, —1 and the upper indices show for which 
operators the (anti) commutator is calculated, while the lower index provides the right choice 
of all Hubbard operators for the r.h.s. of the (anti) commutator. In the case when X- and c- 
subsystems do interact and corresponding orbitals overlap, the constants e in (29,30) contain 
sum of products of matrix elements of the inverse of overlap matrix (determinants) and are 
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given by quite cumbersome expressions. Usually, in physical situation, only some of c- 
electron bands have non-zero overlap with /-electron orbitals because of their location in 
different energy regions. For this reason one needs actually to take into account the overlap 
integrals only with those of /-orbitals which are located in the vicinity of Fermi energy, i. e. 
with upper subset of orbitally polarized /-shell. In order to give feeling of the structure 
of constants e we give a simple example in the Appendix B. These constants can easily be 
calculated for any given set of many-electron states. Below we will assume that ef', etc., 
are calculated and known. One way to find the commutation relations is a straightforward 
calculation by making use of the facts that we know that 

i c jL, fj'mic) = (O'^jLj'miCr, (31) 

and any of the operators X a or h r can be expressed in terms of /, /t— operators. For example, 
using the general expansion 

^ = E(r|^|r')x rr ', (32) 

IT' 

and explicit expression for the operator of annihilation of electron on the third orbital in 
the state with second and third orbitals occupied, we see that 

(1 - ni)^/ 3 (l - n 4 )...(l - n 14 ) = -X^ 23 \ (33) 

This direct way of calculation is quite long. There is, however, a way to make it compact. It 
is easy to verify that one can work with the operator J2jl Oj'njLCjL under an anticommutator 
sign (as one can with /^ (here \i = (mi, a) ). This observation immediately leads to the 
formula 

{cn.,x}}= £ £ (o-V,^/^0 ai 4 1,s 4> ( 34 ) 

aim ; cr £=A,A,F 

which reproduces the one we calculate in a straightforward way from the the expression for 
the X-operator in terms of /, /^-operators. The sign of each term in the sum over a\ is 
determined by the matrix element (f m y) ai , which takes into account the convention about 
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the order of /-operators in a |T) -state automatically. The same argumentation leads to the 
following expression for the commutator with h v \ 

[CjL,h^}= E {O-^LJrn'a'tfrn'v'Y'el^X],. (35) 

In the non-interacting system (Hi nt = 0) inverse processes of excitations are absent since 
{c, /} = 0. This anticommutator, of course, cannot be changed by Hi nt . However, the 
partial anticommutators {c, X a } may be non-zero, if they are considered separately. As it 
will be seen later, these anticommutators may involve two-particle processes. 

Since the majority of rare-earth-based materials have a single-determinant ground state, 
it is quite obvious that the orbital representation is convenient for the description of 
ground-state properties. The excited states, however, are described in most cases by 
many-determinant states (in quantum-chemistry language described as " configurational 
interaction") and, therefore, the procedure of constructing the states |T n ), discussed in 
details by Irkhin and Irkhinll should be used. In these cases the structure constants 
e® b and so on, which come from the commutation relations, and the matrix elements 
{.frrnaY = (fmia) Tn,rn+1 = (T n | f mi a | Tn+i) , include also the Clebsch-Gordon and fractional 
parentage coefficients. Here we will not go here into details of this well-developed theory 
and refer the reader to ref.18. 



B. S-matrix and the exact equations for GFs in functional derivatives. 

1. Quasi- fermions 

Now we introduce the Matsubara's Green's function (GF) for the fermion-like transitions. 
The idea, exploited in Refs.@S0 for the derivation of the diagram technique, is based on the 
possibility to express the terms, which come from interaction in the equation of motion for a 
GF, in terms of functional derivative of GF over external field, and, in the next step, on an 
iteration procedure. The equations of motion for Fermi-like and Bose-like operators are given 
in Appendix C. As seen from these equations (Eqs (C5-C16)), in our case the Hamiltonian 
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generates the terms with the operator Q°; c . Therefore, we have to introduce a type of external 
field U (t) which generates this operator. However, an inspection of the further steps in such 
a version of the theory shows that using the field in the form JUj C (t)Qj C (t)dt introduces 
information in a too much integrated form, which does not provide a possibility to calculate 
the quasi-boson GFs. Instead we introduce the definition for the quasi-fermion GF with 
different fields for different quasi-boson transitions; 

( T Vj\(t)Vj'\'(t )>„ = {T S(-ip,0)) • (36) 

The ^-matrix is then defined as follows: 

S(-i/3,0) = exp{-i r P dt E U x (t)Z*(t) + E c Ut)F jL ,fL'(t)c j/L/ (t)}}. (37) 

a 

The additional terms which are generated by the external field U and which should be added 
into the equations for the GFs are seen from the equation: 

id t {TS(-iP,t)r] jX (t)S(t,0)} = {TS(-iP,t)[idtVjx(t) + A jAj , A ,(t)^ A ,(t)]5(t, 0)}. (38) 

The elements of the matrix A- x -, x ,{t) are given in the Appendix D. 

Second term, c^Fc, has the same operator structure as zero Hamiltonian of conduction 
electrons, c^hc; therefore, we will take them into account if in the equations of motion for 
any operator make exchange h\ 2 — > h u = h 12 + -^12- Actually, we will need this field only 
for consideration of remaining terms of Coulomb interaction. On the other hand, it gives an 
additional convenient freedom to operate with equations for GFs. 

As seen from the equations of motion (see Eqs.C5-C6) for the components of the 77- 
operator, c- and X-operators, the equations for the GFs {Trffj) u contain the Greens function 
(TQ(t)rj(t)fj(t')) u and (Trj(t)Q(t)fj(lf)) u (see Eq.(C10-13)). Both these GFs can be obtained 
from the three-time GF (TQ(t")rf(t)fj(t')) u by taking the corresponding limit, \im t n^ t+0 ... 
or lim t ^ t _ ... . Since Q ah = efZ*, SS/6U^(t) = -iSZ*, and 

(TZ^'HmDu = (<TZ € (O>« + <^;y)<r^(0>„, (39) 

25 



any of the three-time GFs can be expressed via functional derivative of the GF (Tr](t)f)(t')) u 
over corresponding field U^{t"), another trick exploited by SchwingerU and by Baym and 
KadanofiS. After that the equations for the GFs can be written in closed form straight 
forwardly from the equations of motion for c- and AT-operators, given by Eqs.C5 and C6, 
and the relation in Eqn.38. The diagonal expectation values of the diagonal operators, (h T ) 
in the equilibrium state have physical meaning of population numbers of the ion's state V. 
When external time-dependent fields are present, these operators do not commute with the 
Hamiltonian even when the interaction is switched off, therefore, (h r ), strictly speaking, 
cannot be interpreted as population numbers any more. Although it is not completely 
correct, for briefness we will refer loosely to "population numbers" for the expectation values 
of any Bose-like single-site operator. 

For briefness we will not write the complete matrix equation for {r](t)fj) u , instead we will 
specify only the left-hand-side 7/-operator. Then we find that the equation for the set of GFs 
(TcjL(t)f]) u has the following form, 



[hnh^dt ~ (0- l H) 3Wl ]{Tc nLl (t)v)u ~ [A)?J lb (t) + w jLJlb ]{TX b h (t)fj) 



-A 




u 



(40) 



The equation for the set of GFs (TXj 



(t)r))u is 




u 
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x{t%(TXl(t)fj) u + w l jlLl!jc (Tc] lLl (t)fj) u }. (41) 

The equations for (Tc f jlLl (t)fj) u and (TX^(t)fj) u can be obtained from Hermitian conjugation 
of the corresponding equations of motion, 

[H)dt(Tc] iLi mJ n ,S Ll , L - {Tc) iLi {t)n) u {HO- l ) nLwjL ] - (TXl(t)fj) u A^ L (t) 
- (TXl(t)fj)J%(j L (t) = i({c] L ,fj}) u 5(t - t>) + (TXl(t)fj) u w j2 - bJL 

+[<9?(n>.+< s ^ Fy ] 

x{(Tci L ^t)fi) u W hL2jia + {TXlm u tt n }Ul)\0~% M , ]L 

x{W* 2 ,, jlLl (Tc jlLl (t)fj)u + tf 2j ATX^mu}(flnO-% MdL , (42) 

and 

[8 ab {-id t - A ja ) - A*d(TX](t)fj)u - A%(TX»(t)fj) u 
= -i({Xf,fj}) u 6(t-t') 



-l(Q?(t-))u + i 



ba/j.—^ 

3 ■ ,,v 



x{ W j CJiLi (Tc, lLl (t)^) tt + t}*(TX»(0^>„}. (43) 
The products (HO -1 ) and (0 _1 if ) are the matrix products in the cc-subspace only, i.e. 

(HO- 1 )^ = H jWl (0-% LufL , (44) 
The labels ± near the time variable denote limits = t ± + . 

2. Quasi-bosons 



Our recipe for calculation of the GFs is so far not complete; we have to describe how to 
calculate the population numbers and their derivatives. The latter, actually, describe the 
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boson-like excitations in the system. Now we turn to the consideration of this part. 
The "gauge" field coming from the source-term of the S'-matrix is 

A**(t) = et£ut l (t), (45) 

where e^ 1 , similar to Eqn.29, gives the coefficients of expansion [Z^, Z^ 1 ] = e^Z^ 2 . 

Let us write down first the equation for (Zj(t)) u . For that we have to add to the average 
of (54) the contribution from the source term. Then we find 

[8^(id t -At)-Af(t)}(TZl(t)) u = 

E Vik^xATvlM^xMu. (46) 

m Ain2^2 

This equation is valid also for the case of diagonal transition, {Zj~^ T ' T \t)) u = (hj(t)) u . In 
this case Ag = A[r,r] = E r — E r = and at zero fields, U = 0, we come to the case of 
thermodynamical equilibrium , 

niAin2A2 

where the population numbers iVr = (Th?(t)) do not depend on time. Then, the right- 
hand side gives zero identically, since summation in the right-hand side of (47) produces 
the difference of Hermitian adjoints. This equation is useful for checking the approximation 
made but does not give a recipe for finding N r . Therefore, we have to go further and to use 
connection between the GFs 

(TZt(t)2$,(lf)) u = Ktf,(t,lf) (48) 

and their end-factors 

P*(t) = (T[zl(t), Zf,(t)]) u = 5 jr ef,{TZf{t)) u . (49) 

Let us derive equations for the Bose-like GFs. They can be obtained in two ways. One way is 
to use the equations of motion, Eq.(C18) and then to connect the higher-order GF with the 
lower-order one via functional derivative. Another way is to take the functional derivative 
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6/SUf,(t') from both parts of Eq.(46). These ways should lead to equivalent results. The 
equations of motion give: 

[Sf(id t - A € -) - .if (/)7y ;i ,., r! /./') = i6(t - f)Pff!(t) + 

Vlx^iTvlx^VnzxMZf^ ( 5 °) 

Let us consider now the following three-time GF (Tr} n2X2 (t 1 )r)l llXl (t 2 )Zj, (t')) u . On the 
one hand, the two-time GF in the last term of Eq.50 can be expressed as the limit 

t2 li > m o (r^ 2A2 (t)^ iAi (t 2 )4'(^)) u = (Ti^Oi^JO^)). (51) 

of this function. On the other hand, we can again use Eq.39 and express it in terms of the 
GF G = (Tr] n2X2 (t)r]l liXi (t 2 )) u .This gives the following equation for the boson-like GF: 

[(^(OJu + i^l^l^^W^^^Ju, (52) 

where (£°) _1 is defined by the equation: 

{[5f(id t - A f ) - Af>(t)]5(t - t 2 )}L% 2 ^(t 2 ,f) = i5(t - t'). (53) 

The procedure which remains to be discussed is how to calculate the expectation values 
(TZ 3b (t 2 )) u , i.e., population numbers. The natural way to do that is to use the multiplication 
rule for the Hubbard operators X?' 9 ■ Xj' s = \j,p){j,q\ ■ \j,r)(j,s\ = 5 q ^ r X^' s . In the limit 
t' — > t + the GFs K(t, t') or G(t, t') contain a product of two operators at the same moment 
of time and, therefore, we can use the rule of multiplication. In the case of retarded GFs 
an application of the spectral theorem gives an equivalent system of equations in those of 
approximations, in which the analytical extension of the vertex corrections does not involve 
different values of vertexes from different time intervals: 

Nf = [lhn Q (Xt(t)X«(0))} i=j , (54) 

or from Bose-like GFs, 
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Nf = [Km o (Zf(t)zl(0))} i=j , (55) 

where £ = [r,Ti]. 

When the iterative perturbation theory at U(t) ^ is used none of expectation values of 
(Z(t)) is equal to zero. In the limit of thermodynamics (U{t) = ) the non-diagonal numbers 
(Z^) = while the diagonal ones N r and their derivatives ^ r su T ^ u=0 are no ^" ^ e 
derivatives of non-diagonal occupation numbers are charge and magnetic susceptibilities (see 
Eqs 94-95 below). However, the correlators {h r 'h r ) also contribute to both of them. 



3. Comments on the system of equations 

Eqs. 40-43 and 52 are the basic equations of our theory which we will use for generating 
a series of perturbation theory. The iterative terms and lowest approximations will be 
discussed in the next section. Within this formulation of the diagram technique the difference 
between the cases of orthogonal and non-orthogonal bases sets consists only in the number 
of interactions involved and in the definitions of the matrix elements. Contrary to the 
standard case of normal fermionic system the set of equations Eqs. 40-43 for quasi-fermion 
transitions is not yet a complete set even in terms of functional derivatives, as seen from 
Eqs. 41, 52 and will be seen also later. The functional derivatives from expectation values of 
quasi-boson operators introduce also energy-dependent contributions to the self-energy from 
quasi-boson excitations and in addition, broadening of spectral weights. The difference to 
the case of, for instance, standard Fermi-liquid theory is that in the Fermi-liquid the boson 
excitations arise only as collective modes whereas in our case they exist even in the level 
of zero Hamiltonian . At first sight the presence of the GFs (Tc^ 2L2 (t)fj) u for r\ = in 
the system of equations looks somewhat strange since we are considering the normal, not 
superconducting state. Actually, the zero GFs of this kind are indeed equal to zero in the 
normal state; this is simply a manifestation of the fact that our external Bose-like field U 
contains components which describe the transitions [r n ,r n ± 2 ], i.e. transfer of 2 electrons 
simultaneously. The corresponding GF describes process of transformation of two conduction 
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electrons to /-electrons on the same site (or vice versa) and appears only as a perturbational 
contribution. Thus, at zero external field we have (c'c 1 ) = whereas 5(c^ci)/5U may have 
non-zero values. This, obviously, leads to the existence of some non-zero correlations of 
the kind (cic 2 (t)c 3 c\(t')) , and as a consequence superconducting fluctuations above T c due 
to phonons. Some other mechanism which are described by the same function may also 
be influenced by the two-particle intra-atomic correlations (for example, it can make the 
solution with the ci-wave symmetry of the order parameter more favorable compared to the 
one with s-symmetrytiH.) 



IV. APPROXIMATIONS 

Eqs. 40-43 and 52 can be used for iterations straightforwardly. However, since the 
subsystems of conduction electrons and quasi-atomic transitions are coupled, it is more 
convenient in practical calculations to make iterations from a more symmetrical form of 
these equations. Let us first rewrite the equations in terms of the operators T]i,r]\ with the 
indexes 1 = (ji, Ai). In order to do this let us return for a moment to the initial Hamiltonian 
of PAM. All its terms except the one TCq = E*phj can be expressed in terms of operators 

Ka = V2V23V3 = H 23 c\c 3 + W 23 c\X 3 + W; 3 X\c 3 + t 23 X\X 3 . (56) 

The term Til does not cause difficulties since, as seen from Eqns.C5,C6 and C18, in the 
equations of motion it generates a term which also can be written in terms of 77-operators. We 
denote this term as 6\ 2 '■ ^12^2 = -^r^J ]• Therefore, the equations for GFs, Qi2(t,t') = 

(Tr]i(t)r)2(t')) u , can be expressed fully in terms of 77-operators and functional derivatives of 
GFs with respect to external fields, 

(i5 12 d t - A 12 (t) - n° 12 - e r 12 )g 21 ,(t,t') = iv iv (t,t') + 

V 23 (T{ Vl (t),4(t)} V3 (t)4(t')) - (57) 
V 23 {T V i(t 1 ){ Vl (t)Mt)}vUt r )), (58) 
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where Ajlj'l' = 0. Since we have moved the non-/ part of the LDA Hamiltonian from zero 
Hamiltonian to the matrix of interaction V, fl° becomes a purely diagonal matrix: 

^12 = ^l<^12 : ®°ja,j>b — Sjj'SabAa, ^jLJ'L' = ^jj'^LL'^L- (59) 

Here we have included a 5-function into the definition of the time-dependent population 
number, 

VwfaA) = = (60) 

- mwii vUm = s(ti - a^aa). (6i) 

It is seen from Eqns. 40-41 that in order to have a compact form of equations, we have to 
form combinations of the operators {Q)+i-£j with matrix elements, such that they reproduce 
full combinations {771,77!}) {^i,^} i n the equations of motion. The GFs of higher order in 
terms of Eqns. 40 and 41, which are expressed in terms of functional derivatives, can be 
written in compact notations as follows: 

\/ 23 (T{r /l (t),^(t)}r /3 (t)^(t , )) = 

V 23 [V 12 (t + )+R l 12 (t + )}g 31 ,(t,t'), (62) 

V 23 (Tr ] l(t 1 ){r ]l (t 1 ), m (t 1 )}r ]l ,(t[)} = 

V 23 [P 12 (t-) + R[ 2 (t-)}(T4(t)rjl(t')) u . (63) 

The operators R are here defined in such a way that R[ 2 (t + ) = and R{ 2 (t~) = for c- 
electrons, i.e., when 1 = (ji, Li) and 2 = (j 2 , L 2 ). The cases when one of these indexes, 1, 2, 
or both of them describe /-electron transitions, are described with the help of the matrices 
a and (3. They have to provide correct commutation relations and, therefore, we define them 
by the equalities: 

07i, vl} = a&A, { Vl , 7/3} = P?h b Z 3b , {vl 4} = p£i b Z 3b (64) 

where f = A, A, T , ( = [r„,r n+2 ] and ( = [r n ,r n _ 2 ]. The complex index 3 b = (j 3 ,f 3 ) 
in the case of the a-matrix, then, 3& = (j3,(s) i n the case of /^^-matrix , and, at last, 
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3b = (.73, C3) f° r Pil^-mairix. The subindex "6" points out that the index runs over 
Bose-like variables. If, for example, 77 = c, and = X\ we have 

^,k = O u % 3 (f M ref, (65) 

and so on. We will not repeat all these definitions since they are clearly seen from the 
expressions given in Eqs. 40-43. Then the operators R l 12 (t + ) an d -^(O are: 

auo - <4 '™ +0 <,;^) - (66) 

amo . , Hm/i?,,^ , A^ju^y (67) 

where the relations between the matrices a and a, j3 and /3 are as follows; 

Defining in this way a- and /3-matrices in the ^-operators we put in a formal language the 
requirement that R12 = when both ends, 1 and 2, are c-electron ones. Further, as we 
remember, in Eqs. 40-43 (see also here Eqns.78 and 81) the term (Tr]2{r] l ,r]3}(t)r]\,(t')) u is 
expressed in terms of functional derivatives of the GF (Tr^W^i' (£'))«• The matrix of this 
GF consists of the same elements as the matrix (Tr) 2 (t)r)l,(t')) u . Therefore, we can express 
the former GF via the latter one as 

where the matrix r 2 3 makes necessary permutations of the matrix elements and, therefore, 
consists of zeroes and unities only. For example, if we deal with a system which has n 
conduction bands and m intra atomic Fermi-type transitions, the matrix r 2 3 has all matrix 
elements equal to zero except over- and sub-diagonals where all the matrix elements are 
equal to one, 

' r G , 2l,j3«.+m+l) T "(i22,j3n+m+2) ••• ^{j2n+m,j:i2(n+m)) 1, 

T (hn+m+\,j 3 l) = l~(j 2n+m+ 2,j 3 2) = ■■■ = T(j 2 2(n+m),j 3 n+m) = 1- (69) 
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Note that the left (t\ — > tf) and the right ( ti — > ij") limits are included into definitions of 
the operators R[ 2 (t + ) and R{ 2 (t~) correspondently. 

Let us now introduce the auxiliary zero GF by the following equation: 

/ dt 2 [-^5 12 d tl - A 12 (t 1 ) - Q° 2 )<^i - t^vfi^Q = 5 ir 5(t 1 - t[). (70) 

As seen from (59), the matrix fi° 2 has purely local, on-site, characteristics. Therefore, it is 
natural to call the function V 2 y "locator". However, mixing, hopping and overlap connect 
different sites, which leads to dependence of the dressed locators T>(t,t') on the variables of 
different sites; for this reason they will be named "pseudolocators". The inverse matrix of 
V(°\ i.e. [D < i 2 (ti,t2)]~ 1 = X'(7 1 (l,2), is given by the expression in the square brackets on 
the l.h.s of Eq.98. Putting 9 = V = in Eq. (57-58) for the full GF we obtain the equation 
for zero GF Qfii^A)' 

J dt"V^ 2 (tX)Gi>",t') = Vif,(t,t'). (71) 

We will use this equation below for iterations. In this notation the full equation for the GF 
acquires the following form: 

/ dt 3 V^ 13 (t,t 3 )g 31l ((t 3 ,t') - {-9 r 13 + [V 12 (t + ) + R[ 2 (t + )]V 23 
--V24[Pu(t-) + R r u (t-)]r 23 }g 3V (t,t') = VwW). (72) 

The same equation can be presented also in a way close to the standard form of the Dyson 
equation: 

{V^ 12 (t,t 2 ) - S 12 (t,t 2 )}g 21 ,(t 2 ,t') = V u >(t,t). (73) 

Here and below we omit the sign of integration over intermediate time (here t 2 ). Eq. 73 
gives the definition of the self-operator matrix S. However, as seen, this equation does 
not fully coincide with the standard definition of self-energy. The difference between Eq.73 
and the standard Dyson equation consists of the presence of Vw(t,t') in r.h.s., instead of 
Su'5(t, t') and is caused by the non-fermion commutation relations {i]i,r] 2 } ^ 5\ 2 . As will 
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be seen later, these terms give rise to whole generation of new graphs. Part of graphs have 
the same structure that in the theory of fermions. The S'-self-operator contains both type of 
graphs, of Fermi-liquid type, and the terms, generated by the right-hand side of Eqs. 72-73. 
In order to distinguish between them we use the term "self-operator" in definition for the 
full GF(23), while the term "self-energy" will be used in a more traditional way, for the set 
of graphs which renormalize the denominators of the pseudolocators, and will be introduced 
later. In order to separate these two type of graphs we will search the solutions of (73) in 
the form (71), i.e. we introduce the full pseudolocator V by the Dyson equation 

{V } 2 (t,t 2 ) - S 12 (t,t 2 )}V 21 ,(t 2 ,t') = 5 u >S(t-t'). (74) 

Then the full GF is given by 

Sir(M') = V 12 (t,t")V 2V (t",t') = V 12 (t,t')V 2 y(t'). (75) 

(Note, that the population number V in the equation above depends on the right time t'). 

The lowest approximation as well as the equation for the self-operator, which is con- 
venient for iterations, can be obtained from its explicit expression in terms of functional 
derivatives. The equation for the self-energy follows from a comparison of the two forms of 
equations for the GF, one, Eqn.73, in terms of self-energy, and another, in terms of functional 
derivatives Eqn.72, and we obtain, 

V(M") = -0 T n// S(t - t") + {-[P 12 (t+) + R[ 2 (t+)]V 23 

--v 24 [P u (r) + R^iin^ygsrit^yg^iit'X). (76) 

A. The Hubbard-I approximation. 

The Hubbard-I approximation is obtained if we put all the functional derivatives equal 
to zero, RQ = 0, which results in 

S$(tX) = = " + V ^ t+ ) V ^" ~ V24pu(t-)T 21 „}. (77) 
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Since this term does not contain the factor the S'-self-operator here can be interpreted 
also as self-energy (denoted as £). If this expression will be used for a system in thermo- 
dynamical equilibrium and far from a phase transition to a state with intermediate valence, 
Hubbard's arguments can be applied, namely, that instead of fully self consistent calculation 
of population numbers one can use just the atomic-like solution for them, Vuitf ) « Vvt(tt)- 
Note that by putting functional derivatives of the GF equal to zero, we automatically neglect 
all the graphs, which renormalize the population numbers. Thus, this is the lowest possible 
approximation. In the equilibrium state external fields U = 0. Therefore, all non-diagonal 
population numbers Vy[ =0. Then the corners (Tec), (TXX) of the full matrix of GF give 
zeroes and it is sufficient to consider only the sub-block with i] = c,X and = c^,X^ . 
Then the expression 



4?(*i,* 3 ) |«=o= -s(tt -hM 3 + v 12 (tt)v 23 ], 



(78) 



gives the following system of equations in the supermatrix form, 

^ ADcdtyjtnL! - ®jL,j 2 cW* 2CdlLl -w jLJlb - 9 jL ,j 2C t cb 



3,31 



V 



-w 



ja,jiLi 



m x (t)]ji ilb -i! ab 



X 



3,31 ) 



V 



Gjib,j'L'(t Gj lb jiai(t t') 



J 



( 



iS(t - 1') 



o- 1 

U jL,j'L> 



A 



r>abr\ _ r ryaa' 

\ ' j u jb,j'L> <\i.ri y ■ ) 



(79) 



The Latin letters D are used for cc- and XX-components of the pseudolocators, whereas 
their inverse pseudolocators here are the following differential operators: 



*lDcc(t)]jl, fL , = 5 jf 5 LL ,id t - (O"^, 
i[D xx (t)]-l fb = 5 jf 5 ab (id t -A,), 



(80) 
(81) 



and the interactions which are renormalized by (convoluted with) the population numbers, 
are denoted by a tilde, 
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W* = P cb W*- i cb = P ca f ab 7 T,T — pab T f\ ab — pacflcb /oo\ 

Finally, the G-matrix is: 

©jLJ'c = O j lj, M (f M ) C , ®jb,j'L> = (fM) b OjM,j>L>- (83) 

Thus, the new element, which appears in our Eqs.82 compared to normal LDA-based cal- 
culation, is that both effective mixing and hopping matrix elements are decreased by the 
population numbers of many-electron states of the ion, since usually they are less than one. 
The latter is always the case if mixing of localized and delocalized electrons is not equal to 
zero. It is worth noticing that any further corrections to the population numbers, which 
do not change their locality and do not introduce an energy dependence, do not change 
the structure of these equations. The approximation, where these population numbers are 
dressed, but still do not depend on energy, we will call "Hubbard-I" approximation. The 
corresponding self-energy is shown in Fig.l. 

Let us emphasize what is achieved for a moment by this construction, i.e., which are the 
new features of our " Hubbard-I" equations, which differ them from the standard LDA ones. 
First, we passed from the ground-state description of the system to the statistical one, 
therefore, some of thermodynamical properties can be calculated. 

Second, we know now where the /-level is situated with respect to the Fermi energy and 
other bands just by construction of the Hamiltonian. This is important for example, for 
calculation of contribution of mixing interaction into effective interaction of conduction and 
localized electrons. 

Third, in spite of the allowance of mixing interaction between localized and delocalized 
electrons the atomic-like Hund's rules can be taken into account (already at the level of 
constructing of the zero Hamiltonian ). 

Fourth, the possible degeneracy of many-electron states involved is taken into account by 
the form (25) of the population numbers. 

Fifth, since the single-site contributions from the Coulomb and exchange interactions are 
included by construction into the energies of many-electron states Er, the calculation of 
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an effective exchange interaction between the angular momenta (also, interaction between 
quadrupole momenta, etc.,) of different atoms does not require special care, which is 
necessary within the standard LDA model for heavy rare earthsS 

Finally, the interaction of quasilocalized electrons with itselves (self-interaction correction) 
is taken into account by construction. 

It is obvious from a physical point of view, however, that the approximation described 
above has very serious drawbacks. If some of /-ions has n /-electrons in the ground state, 
the state with n + 1 electrons belongs to the excited ones and usually is separated by such 
a large energy gap from the ground state, that the energy difference Ap"^' n = E^ 1 ^ — Ejfi 
exceeds the Fermi energy and lies in the energy interval of the delocalized part of the spec- 
trum of single-electron excitations. Corresponding excitations unavoidably should become 
delocalized. Therefore, part of the spectral weight should be moved into this energy region 
and the renormalization of the population numbers (combination of which describes in this 
approximation also the spectral weights) and, also, the energies of the transitions Ap^' n , 
should contain the contribution from this delocalized part of spectrum. The Hubbard-I 
approximation misses it. These corrections arise only in higher approximations. The other 
drawback of this approximation is that the atomic-like /-levels remain non-renormalized 
by the interaction of the /-ion with surrounding medium. This shift should arise due to 
admixture of the wave functions of neighboring sites to the one of the reference site and 
interactions of /-electrons with bosonic excitations (analogously to the Lamb shift in quan- 
tum electrodynamics). All these effects can be taken into account via the corrections to the 
self-energy and we turn now to the formulation of this part of theory. 

B. Vertexes and equation for self-operator in terms of vertex. 

In order to obtain the next corrections we insert the GF, Q 36 , in the form of Eqn.75 into 
the expression for the self-operator, Eqn.76. The differential operators R acting on the GF 
Q gives, 
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iWti)£3i'(M') = Rii{h)[v^{tJ)Vsv{t')\ = 
[Rnit^V^t'^Vwit') +V 36 {t,t')[R 12 {t 1 )V 6 v{t')} = 

-v 39 {t,t 9 )[R 12 {t 1 )v 98 \t 9 ,t 8 )}g 81 ,{t 8 ,t') + 

T> 3B (t,l/)[R Vi (t l )V 61t (l/)] (84) 



where the trie" 



Ri2(ti)[V^{t 3 ,t 6 )V 67 (t 6 ,t 7 )} = => 
RMV^te) = -V^tg^R^V^itg^V^te), (85) 

has been used. Having in mind that the i?-operator is equal to zero when both indexes run 
over c-electrons, we restore its definition via the matrices a for the R- and f3 for the R r - 
operators and introduce the vertexes in a standard way, namely, as functional derivative of 
the inverse GF with respect to external field, 

R l Att) Vnfah) = f lim a%4 V 98 l (t 9 ,t 8 ) 

t4-»ti+u vMfa (m) 

lim 0(12,4^984 (tg,t S ;U) = 0(i2aT 1 7S a (t 7 , t$] tf), (86) 



Ru(tl) V,i(t 7 , t 8 ) = lim J^, 3b — V£(t 9 , t s 



= t J*^_ Q ^12,4^98,4;, (*9, ^ h) = (3^ Ab T r 98Ab (t 9 , t 8 , t x ). (87) 

It is worth noting that in our case the role of standard GF is played by the pseudolocators, 
and in this respect the introduced vertexes are indeed the standard ones. The new elements, 
which arise in the techniques for non-fermion operators, are the derivatives RV . They can 
be written via the correlation functions of Bose-like transitions (see Eq.95 below) and we 
will discuss them a little later. The other new element (compared to the standard case) is 
that the vertexes which enter the equations in the " left" and "right" limits. The origin of 
this difference is seen from the definitions of the operators Ry^X') an d ^12(^1"); the vertex 
in the left limit arises in the usual technique too, while the one in the right limit arises from 
transitions which involve creation or destruction of two particles on one site simultaneously. 
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This process is described in the equations of motion by the anticommutator {r/i, r] 2 }, which is 
equal to zero in the case of usual fermions, but is not in the case of two-electron intra-atomic 
transitions. 

Using these expressions we obtain the following equation for the self-operator, 

Sn»(M") = - W23V 39 (t,t 9 )a% b T 9VI , 5b (t 9 ,t'';t + ) + 

-V 24: T 23 V 39 (t, t 9 )(3i2 bb T 9V , bb (t 9 , t"; t ) + 
(-V2sV 36 (t,t')[R[ 2 (t + )V 6 At')}- 

-V 24 T 23 V m (t,t')[R r u (t-)Vei>(t')} ) Gy\„{t',t"). (88) 
% 

Thus, structurally the self-operator can be separated into the sum of a few terms 

S = S H + S 1 + S P , Sh = E#; Si = Ei; E 1 = E L + E^, (89) 

where E L and E^ are the terms which contain left (~ a) and right ((3) vertexes and are 
genius self-energies E (i.e., the graphs, which cannot be separated into two non-linked 
graphs but cutting one line of D-function), while S P , the contribution, proportional to 
G^ 1 , contains derivatives of the population numbers and renormalizes end-factors. Let us 
also consider the corrections to E L First we need the bare vertex. The definition of it arises 
if we replace the full inverse pseudolocator in Eqns. 86-87 by the bare one, V^ 1 — > Vq 1 . Since 
R^^Vq 1 = a(5T>Q l / 5U) = —a(5A/5U), and a similar relation holds for R r (t)V, the number 
of different bare vertexes is determined by the number of different combinations of external 
fields A, (Appendix D, Eqs.1-8). From the definition of the self-operator V' 1 = Vq 1 — S we 
see that the exact total vertex being derivatives of the self-operator with respect to external 
fields, can be rewritten in the form of a sum of the bare vertex 7 and corrections, 

^98,56^9^8^5) = 798,5 6 (*9>*8; h) , (90) 

^s,^ts-,t 5 ) = -5(t 9 -t 8 ) 5 -^^. (91) 
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These relations are useful for obtaining subsequent corrections to the vertex, substituting 
into Eqs.86,87 approximate expressions instead of the full self-operator. We will return to 
it a bit later. Before that, let us rewrite with the help of Eqs. 76,88 the equation for the 
full GF in a form which clarifies the physical meaning of different terms in the self-operator 
(Eqn.89). Using the pseudolocator in the Hubbard-I approximation, V^ 1 = £>q 1 — H H , we 
find, 

{(P"(M"))n" - ip 3 9(t,t9)[U 2 3< + 5t r 91 ,, 5ft (t 9 ,t";t + ) - 

- Vmtw^iv 5b (t 9 , t"; t-)]}g v , v (t" , t') = 

P ir (M') + APir(M'), (92) 
where the renormalization of the population numbers is given by the following expression: 

AV lv (t,t') = --y 23 V 36 (t,t')[R l 12 (t + )V ev (t')] + 

-V24T23T> 36 (t,t')[R[ 4 (t-)V 6 r(t% (93) 
% 

Thus, we see that the self-operator can be represented as S = £ + AP ■ G^ 1 what allows 
to separate the graphs to the ones which renormalize the pseudolocator D and are the 
contributions to the self-energy, and the ones which renormalize the end-factor AP. Let us 
turn now to the Bose-type of GFs. In order to obtain the equation for these GFs, we will 
use again the short-index notations 1& = (ji£i) and rewrite the equation for the boson-like 
GF in terms of vertex, Eqn.86 (we have to use once more time the Kadanoff and Baym's 
trick-Eqn.85). Then we obtain, 

[L'ihM^K^M) = S(h - O^WO + 
hm J 7 V 1 1 £(TZ 1 ,(t')) u G 21 (t u t 2 ) + 

t lim n ^[- J D 24 (ti,t4)r 4 5;i'(t4,t5;t , 1 )G' 51 (t 5 ,t 2 ) - 

D^tut^a^K^^U) ~ (TZ 1/b (t[)) u (TZ 3b (t 2 )) u ]. (94) 

Here we have used the definition of Psvfo), Eqs. 60 and 61 and represented the derivative 
SP/5U in the last term of this equation in the form, 
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= a 
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{ 3 MK 3b;1 ,(t[,t 2 ) - (TZ 1/b (t[)) u (TZ 3b (t 2 )) u }. 

t 



(95) 



This completes the derivation of the exact equation for the boson-like GFs. Eqs.94,95 
should be used also in Eqs. 84,88 and 90 for calculation of R l P and R r P. 

The way of deriving different approximations for the boson type of GFs is quite similar 
to the one used above for the fermion-like excitations. First, of course, we have to define 
the zero-order GF, putting V = in Eq.94, which gives 



Then an analogue of the Hubbard-I approximation for fermion-like transitions arises if we 
put in Eq.94 the vertex T = 0. As seen from the Hamiltonian, within the Hubbard-Anderson 
model, the only way a renormalization of the boson-like GFs arises is due to mixing inter- 
action and hopping, therefore, it requires knowledge of the electron GF. The latter creates 
fields which shift the energies of the bosonic transitions. 

C. Mean field approximation and self-consistent equation for Hubbard U. 

As seen from Eq.92, the part of the self-consistency which does not involve the effects of 
Sp, can be provided in the usual way: a choice of certain approximation for the vertexes, 
r — > Y^ appr \ leads to a corresponding expression for the self-energy via Eq.(88), while the 
GFs entering the Y^ appr ^ (if r ^ 7) contain the same £ and r( appr ). According to this scheme 
the simplest approximation which is beyond the Hubbard- 1 approximation, should contain 
the vertex. It arises if we neglect the terms RT, in the equations for T, i.e., r — > 7. Since the 
//-part of Dq 1 is purely local, and the bare vertex does not depend on energy, the diagonal 
parts of the loops 
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£i7(M") = -V 2 3^39(t,t9X + 5 b 7 91 »,5,(^,t";t + ) + 

% 

-v 24 r 23 v 39 (t, t 9 )/3EU79i'',5 6 (*9, n 



(97) 



renormalize the centers of /-bands. 
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Indeed, since we assumed that /-functions, belonging to different sites do not overlap 
with each other, O i) _ L= /y i i/=/ = 0, we see that the energy A s of the transition a = [T n , T n+1 ] 
is shifted by 5T} oop ; if we neglect the non-diagonal corrections, we find the following equation: 

A* = A, + (E^) 1=M) , 2=M) = A a + SA- a {Al). (98) 

This is self-consistent equation for transition energy; renormalization comes from hopping 
and mixing interaction. If we take Hubbard's estimate for E^ n , E^ n = e°j-n + ^Un(n — 1), 
then the energies of the transitions A° = Aa=[r n+1 ,r n ] = E rn+1 — E Tn ~ U/2; the lower 
transition A® = A a= [r n ,r n _i] ~ U/2 too. Here we used the Hubbard's estimated! for 6®: 
the main contribution to this energy comes from attraction to nucleous, which provides 
localization of certain number no of /-electrons. Therefore, dE n /dn = gives e° = — (uq — 
\)U. Correspondently, a fermion //-GF has two poles. The lower one is located near energy 
of lower transition A° and the upper one A° ; the Hubbard bands are developed near these 
energies. The shifts of A° (i=l,2) due to hopping and mixing can be roughly estimated as 

V*" i (n)V ai (<n') 

^A^EipW^ + E \ 1 7 {P> f(E;)^ = l,2. (99) 

Here t l hop {p) is an effective hopping, containing all chain (~ Hubbard-I approximation) con- 
tributions. However, at large enough A main contribution comes from the second order 
(~t 2 /A~t 2 /f/). 

It is also easy to see that the equation A* = A + <5A(A*) can have two solutions, one of 
which is located near Fermi energy, which is often called as "Kondo resonance"!^. These two 
equations show that, contrary to the s-band Hubbard model, one can expect that in reality 
the renormalizations of A\ and A° should be quite different, at least, for rare earth ions: on 
the one hand, hopping and mixing are strongly suppressed in the region of energies uj ~ A\ 
since this energy is very deep under Fermi surface and even below the bottoms of conduction 
bands; on the other hand, in the region of energies to ~ A° hopping and mixing should be 
much stronger, therefore, in the region of energies in vicinity of Fermi level the mechanism 
described by Eq.99 can be much more effective. It is clear, however, that this mechanism 
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is not able to provide strong reduction of the bare Hubbard repulsion A° — ~ U from 
20eV to 5 -T- 7eV (the latter is the value usually obtained within LDA-based calculations). 
Therefore, other, stronger contributions, should be generated by full Hamiltonian. Never- 
theless, it worth to note that the non-diagonal matrix elements of hopping, t XfJ "(f{) a (f\) b 
and effective hopping via mixing, ~ [V^ ai (p)V bi (p) / (00 — £\ p )](l — <5 a b) , split the Eq.99 to the 
system of equations. 

At last, its c/-parts renormalize the mixing interaction. The self-energy which includes 
the Hubbard-I and the loop correction Eqs. (97,98), still does not depend on energy. Since 
this approximation is also beyond the Hubbard-I approximation, we will call it mean field 
approximation (MFA). Similar approximation for multiple-orbital Hubbard model but for 
the simpler case of orthogonal basis set has been considered earlieiHl. The fact, that the 
equations for fermion-like GFs always contain vertexes in combination with the matrices a 
and j3 allows to use a shorter graphical notation: the sums a ■ 7 and (3 ■ 7 can be represented 
as four-leg vertexes, as shown in Fig.2. Then, in this graphical language the loop correction, 
Eq.97, acquires the form shown in Fig.3. 

D. Further corrections to vertexes and self operator. 

In order to obtain the second-order corrections to vertexes and self-energy we insert = 
5Y> loop , Eq.97, into the definition of the vertexes Eqs. 86, 87,90, 91. Since neither the interaction 
V, nor the matrix r depend on the external fields, only pseudolocators entering the E'- 1 *' give 
non-zero contribution when we calculate the functional derivatives RT,. Therefore, we can 
construct the higher-order corrections simply by formulation of the rule of differentiation of 
the pseudolocator lines. This rule is displayed in Fig.4. Applying it to the graphs of the T.^ 
, shown in Fig.3, we obtain the desired corrections of second order to vertexes (displayed 
in Fig. 5). Then, inserting these corrections to the definition of the self-operator via the 
vertex, Eq.88, we obtain the corresponding second-order corrections to the self-operator. 
The graphs for these corrections are shown in Fig.6. The corrections to the vertex in Fig. 5 
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describe contributions of electron-hole excitations to the kinematic interactions , which arises 
due to strong correlations and are caused by mixing and hopping. Note, that there are also 
other, formal, difference with the standard case: instead of pure GF one of the shoulders of 
the electron-hole loop contains the product of the line of interaction V (wiggle line) and the 
pseudolocator. 



E. The screening 

The next obvious step is to take into account the screening effects. The lowest approxi- 
mation is an analogue of the standard random-phase approximation, used in the theory of 
normal metals. In order to avoid possible confusion, we remind the reader that here again 
it renormalizes not the Coulomb interaction but only mixing and //-hopping. One way 
to construct it would be to solve the integral equation for vertex, generating the series of 
"electron-hole" loops (see Fig.4) and insert the solution of this equation, i.e., corrected ver- 
tex, into self-energy for pseudolocator; then to derive similar corrections to the population 
numbers (r.h.s. of Eq.92 for the GF). The equivalent, but shorter way, which automatically 
renormalizes the interaction V in all corrections, is to make an exact transformation of the 
shift in the external field of Eqs. 40-43,52,92,94, introducing the description in terms of ef- 
fective fields, $ = A + YS loop \ The quasifermion GFs depend on the external field U only via 
the "gauge" field A 12 (t) = v^U^t) , where the matrices v are listed in Appendix D. There- 
fore, we can express equations for these GFs in terms of derivatives S/SA e ^, schematically, 
we have 

SU h (t) [ SU h (t) ; 4 12 (fi) SUu(t) *i2<«) 

Let us now perform this calculation. First we have to separate out the loop contribution (97) 
of the self-energy in the equation for GF (Eqn.92). This one- loop correction contains zero 
vertex 7; self-consistency requires that the full GF, which enters the loop, should contain the 
same self-energy, which is taken into account in the equation for GF. Taking into account 
that 
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79i»,5 6 (f9,*";f) = -S(t 9 - t")5(t 9 - t)u 5 9 b r/ , (101) 
and including the self-energy correction S^ oop ) into the effective field, 

$n»(M") = A n „(t,t") + -a^ 5b V 23 D 39 (tt + y 9 b vr 5(t + - t") + 

-P?2zy2iT 2Z Dw{tt-)4{,,5{t- - 1"), (102) 



we find, that the factor Yiyi 2b (t,t') = 5§nr/(t)/5U 2b (t') which renormalizes the interaction, 
should satisfy the following equation, 

Y lv .,2 b (t,t') = v 5 9 \„[5 91 5 2b5b 5(t - - D ^Mt,t 10 ) S -^^y^ x 
{a^ 5b V 23 D 11 , 9 (t 11 ,t + )5(t + -t") + 

PT4- 5b V24T2sD 11;9 (t 11 ,t-)5(t- - t")}\. (103) 
Eq.103 contains full pseudolocator D. 

1. Fermi-liquid-like random phase approximation (RPA) 

If we substitute into Eq.103 the bare vertex 7 = SDq 1 jSU\ h instead of the full vertex 
T = 5D~ l /8U\ b , we obtain only the first correction to Y. In order to obtain a closed equation 
for Y, we have to take into account in D the same correction to E, namely, Tf oop . Since 

5D^/5U lb = -5<$>/5U lb = —Y, (104) 

with £)q/ = D^j 1 — E^ oop ) (5-functions are omitted) we obtain the desired equation, 

yn»,2 6 (M') = is 5 9 b A5 9 i5 2b5b 5(t-t') - D 3tl0 (t,t w )Y w , 11;2b (t 10 ,t') x 
{a^ 5b V 23 D 11 , 9 (t 10 ,t + )S(t + -t") + 

^5^24X23^11,9^10, i")5(r - if')}]. (105) 

Note that the complex structure of coupling of the Fermi-like excitations to the Bose-like 
ones makes the screening matrix dependent on three indices. 
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2. Full equation for the GF in terms of RPA GF 



Now let us write down the self-energy part of the self-operator in the l.h.s. of Eq.112, 
which can be named Ed (the index D here emphasizes that it relates only to the pseu- 
dolocator D), in the form: = H H + Y} oop + Y7 est . The effective field in these terms is 
$ = A + Y} oop . Note that since is not local, it cannot be included into effective field <3>. 
Then, the vertex T, defined by Eqs. (86), (87) and (90) can be expressed in terms of Y as 
follows, 

r 9 i'', 56 (*9,*"; u) = -[Y 9V >, 5b (t 9 ,u)5(t" - * 9 ) + 91 "\ 9 \ J y 67i5b (t 6 ,t 4 )]. (106) 

With the help of this expression for the vertex we can represent the equation for the full GF 
(92) in a form which does not contain the non-screened interaction, 

P*kn"(M") " ±D Sfi (t, t 9 )[a^ 5b } V 23 Y 67 , 5b (t 6 , t + ) - 

Xyrest(f f//\ 

flZ^ J f^g W 2 ^ b {t^r)\}G vl , v {t'\t') = 

P u >(l,l') + AP n ,(l,l'), (107) 

where the derivatives 5/SU in the formula (93) for AP are also rewritten in the form 5/ 514 
= Y 8/8$. At last, the equation for the pseudo locator D RPA becomes, 

[5{t,t")5 12 id t - Ef 2 (M") + -V 23 a^ 5b D^ 9 FA (tX)Y 9 r', 5b (h,t + ) - 
-V24T 23 (3TA, 5b D^ FA (t, t")Y 91 „, 5b (t 9 , t-)]D RPAM (t", = S 2 ,5(t" - t'). (108) 

The role played by Y in the strong-coupling theory is analogous to the one played by dielec- 
tric function in the weak-coupling theory of Coulomb systems: it describes the screening of 
the interaction V 12 . 

V. COULOMB INTERACTION 

The detailed consideration of the Coulomb interaction is given in next paper, here show 
the technical idea of description and emphasize the new elements arising within this ap- 
proach. 
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Full Coulomb interaction have been written in X-representation as follows: 

h cou1 = £W4 + (fl) ai x?)(4 + {flf*x?){c z + (f 3 rx?)(c 4 + (h) aA xr). (109) 

The single-site terms vppff have been transformed into X^-Er^J; the terms which 
contain three single-site /-operator and one operator of conduction electron are included to 
mixing interaction Vc^X + h.c; the terms of the kind v(fi p f)jf /(l — ^y) are included into 
Hubbard's hopping Hamiltonian. Note that in case if some selection rules happen to be valid 
for the single-particle matrix elements of hopping and mixing interaction, the effective matrix 
elements have different selection rules due to presence of Coulomb terms in it. All the other 
terms should be considered. Remind that due to large strength of on-site interaction we are 
not allowed to decouple the operators belonging to the same site. Thus, if these /-operators 
are not neighbors in the Hamiltonian, we have to move (using rules for anti commutations) 
them to each other and then either to transform the /-product into Hubbard operator(s) 
or to transform each /-operator into X-operator and multiply them. This procedure also 
helps to clarify a physical meaning of corresponding terms. The number of different terms in 
the Hamiltonian is large and we are not able to give here all formulas. Besides, a degree of 
importance of different terms is material- and scenario- dependent. Thus, we will consider 
here few terms in order to show how the developed above perturbation theory is to be 
modified in order to include Coulomb interaction between different sites (the description of 
the Coulomb intersite interaction within slave-boson technique has been considered earlier 
in Ref.3). 

Let us start with the terms, which describe simultaneous "hopping" of two electrons 
from one site to another: 

h? = E^mMfif^(fshfzj4 = £cf 44- (no) 

One can expect that this term can be neglected in the lowest approximations since it involves 
large-energy transition (Hubbard U in denominator): although [X a , Z v ] gives X b in equation 
of motion remaining operator Z v requires in order to obtain non-zero result; the GF 
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(TZ v Z f ') gives makes the contribution small. The reason why these processes have to be 
considered is that there are high-order processes (like assisted spin-flips) which contain in 
denominator small difference of these large energies. 
The terms 



UlJJ{fljJ\zl4 = £c$44- (in) 



do not change number of /-electrons in ions. Part of these terms describe Coulomb screening 
of ions, part can be transformed into spin-spin, orbital-orbital interactionsill and so on. We 
will not consider particular terms from the sum over £, £' here, however, note that on one 
hand, instead of one Heisenberg term — J^/JjJj/ the expansion 

H exch = -^/^ 7 (J 1 R) Q (J 2 R) /3 (J 1 J 2 )VR a+/3 (112) 

arises^ (with a + 7 < 2l± + 1, (3 + 7 < 2l 2 + 1, a + f3 being odd and vector R connecting 
ions); on the other hand, the different-site-exchange integrals Vj^f^j^jr^ are small due to 
small overlap of /-functions belonging to different sites, what makes the exchange-type of 
contributions from mixing and hopping visible on the background of positive direct exchange 
integrals and make non- ferromagnetic ordering (e.g., antiferromagnetism) possible. Now let 
us write down the terms, which are added into the equations of motion for the functions 
(Tc n irf) u : 

Jim 0£ Jl fiE!*C$f[(zf l {1»)) +i_^- T ]<TA i W(f)) u + 

r 
8 



i"-< ' '■' " ' dUf,(t" 



^"i^ (113) 

They contribute to the self-operator Tj c xGx v - The structure of this correction is similar to 
the term which comes in Eq.78 from hopping. Obviously, the Hartree corrections which is 
obtained by putting 5(Tcq^) u /5U = 0, are the biggest terms, and of these terms the biggest 
one is the diagonal term, with £ = [r, T]. Next terms, generated by 5(Tcrf) u /5U, in lowest 
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approximation give the loop-corrections to the dielectric function and depend on the scenario 
under study. 

Let us obtain now the contribution from this interaction to the equation for (TX a r) ji ) u . 
Calculating corresponding commutator, we find quite similar result, namely, the term which 
should be added to T^xxGxrj, has the form: 



t lim e?Cp nj [(Z?{t")) + i-^r^WCf )>« + 

This term has the structure similar to the one ~ Q£ in Eq.40. In short notations we can 
write the two ones above as 



^ + ^ 2 ^i b2 J(^iJt / 0)+^^^](^ 2 (t)r / t(t')) u , (115) 

where notations are obvious from comparison with explicit form. Thus, these terms of Hamil- 
tonian do not generate essentially new elements in the digram technique. The Hartree term 
eC(Z} u (the terms at eC[ 5/5U] (TX a rf) u = 0) describes the static screening of charges of 
neighboring ions, (Z — p core ) — p f . The non-spherical part of it is responsible also for the crys- 
tal field splitting of the states Usually this splitting is less important for the Fermi-like 
excitations than for Bose-like, where it can be seen directly in neutron scattering experi- 
ments. Sometimes, however, a scattering of conduction electrons on Bose-like excitations 
in crystal electric field can strongly renormalize effective electronic masses@0. The terms 
~ eC[ 5/5U](TX a rf) u describe effects of exchange and the dynamical screening from those 
/-bands which cross Fermi energy. The /-density in the close vicinity of Fermi level Bp in 
normal materials is usually small and, therefore, these terms can be neglected. These contri- 
butions may play an important role in physics of the intermediate- valence (IV) compounds^ 
and, also, of rare earth metals0. In both latter cases the /-level can be situated just slightly 
above Ep. It can be important also in d — (^-electron scattering in Fe,Co, Ni metals and 
their compounds. 
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Let us pass to a next term, v 1234 a\a\ a^a 4 . Again we have to calculate the corrections 
which arise in equations of motion for GFs. Of course, all terms with two /-operators and 
two c-operators should be taken into account. Since they give similar contributions, for 
briefness we consider only one of them, 

VnLJuju'n' L'{^ n L c n> L'X/j/J ' Z j = C nLn> L'i C nL c n> L>)Zj . (116) 

Summation over transitions £ includes both, the Coulomb interaction of the density- 
density type, and exchange interactions. Commutation of c-operator with this term gen- 
erates in the equation for the GFs (TciL(t)rf(t')) u the GFs: (Tc n2 L 2 (t)Zj(t)rf (t')) u and 
{Tc^ niLl (t)c n2L2 (t)X b j(t)rf(t')) u . One can express them via functional derivatives of GFs as 
follows: 

}}^o-% L , niLl ci\ Lin2L ^ 



3 

^n 2 L 2 (*)>« + 



\ F 6 ( ,J (TX](tW{f)) u . (117) 

""n\L\n 2 L 2 \~ ) 

First line gives contribution to S cc G crj ;the Hartree term, ~ 0~ l C {T Z) , describes a simple 
fact, that the conduction electron "sees" screened by /-electrons ion's charge, (Z—p core )—pf, 
as well as with average angular moment if the phase under consideration has a long-range 
magnetic order. At zero external fields (TZ) — > (h r ) = iV r and (TJj) — > (JJ). Second 
line in this formula contains new element, S/SF, describing creation of the conduction-band 
electron- hole pair (contribution to S c xGxr))- 

Let us turn now to the GFs {T X^{t)r^ {t' )) u . From equations of motion we find that the 
correction to right-hand side, 

Cn lLl n 2 L 2 X ($ G Iv)nlL 1 n 2 L2 ( 118 ) 

involves quite complex GFs: 

(SG Iv y n % in2L2 = (0-X,n lLl (flM(Tzf (t)c n2L2 (t)Zl(t)^(t')) u - 
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(O- 1 )n a L a ,^«Tci l£l (0^(0^(0»7 t (O>« + 

The GFs in first two lines contain per two Z-operators, therefore, in our equations appears 
the new element, derivatives of second order: 

{TZf (t)c n2L2 (t)Zl(t)r)Ht')) u = 

t hm + [(TZf(t")) u + i^^](Tc n2L2 (t)Zl(t)rj\t')) u = 

lim [(TZf(*"))u + i— J ] x 

{lim [<r^(f 3 )>„ + i— |— KTc^f^WU- (120) 
Thus, we've obtained the correction to £xcG CT; . In the same fashion we find that 
(Tcl Ll (t)ZUt)Zl(trf(t>)) u = 



(£ n[(T2?(0>. + i 5 ^]x 

{lim |(TZj(( 3 )}„ + i *-^]<Tc* lil (t),'(f )>.}• (121) 

Note, that this term does not give corrections of mean-field type, of first order, ~ v J n \ Lin2Li . 
At last, the function (T(JcX)rf) u can be written down as 



(^(0^(^(0^(0).= 

lim[(Tct iLi (^)c n2L2 (^))« + - , H M (t)JV))u 

and gives contribution to T, X xGx v - 

The mean-field contribution from the last term, 5A^ oul = n2 L 2 (7 1 c^ lLl c n2 L 2 ) u £^, 
should be added to the equation for effective energies of transitions A^ and Ag (see 
Eqs.98,99): 

A* = A° + A™ 9 + A% oppin9 + Af ou, J 2 = 1, 2. (122) 
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As we discussed above, the correction /\™ lxina _|_ popping ^ Q ^ e deep level A^ are small due 
to strong suppression of mixing and hopping in this region of energy. The corrections from 
vc^cZ are not so sensitive to the position of the A°-level and depends mainly in the density 
of c-charge or, since the contribution is local, just on number of c-electrons per /-cite. In 
order to visualize (see discussion in Ref.0) how this correction works, let us write it for a 
moment in the form Vn c rif. Then E^ n = £o n / + ^ n f( n f ~ -0 + Vnjif and we find that 
Ap n+1 r n = el + Unf + Vn c , A rra r n _ 1 = e{ + U{nf — 1) + Vn c . Thus, the positions of both levels 
are essentially shifted by the Coulomb interaction with conduction electrons towards Fermi 
energy. But the effective Hubbard U* remains unchanged, U* = Ar n+1 r n — Ar n r n _i = U, 
until number of condition electrons n c remains unchanged. Therefore, main mechanism of 
a decrease of Hubbard U is expected from the processes of dynamical screening or, in other 
words, from partial derealization of wave functions of /-electrons. The calculations, based 
on the density-functional theory take into account these processes very efficiently (even 
slightly overestimate) giving, thus, much smaller values of U* than the bare value of U. 
Accurate evaluation of screening effects requires special analysis and we postpone it for a 
future work. 



VI. DISCUSSION AND CONCLUSIONS 

The theory presented here is intended mainly for those systems which contain a group of 
correlated electrons in the strong- coupling regime. The correlations have local, intra-atomic 
nature and therefore it is reasonable to construct a perturbation theory for the description 
of these electrons starting from the atomic limit. The technical difficulties in our approach 
come from the necessity to perform renormalization of fermion-fermion interactions in the 
site representation, which unavoidably leads to non-fermionic commutation relations of new 
on-site collective variables. This, in turn, makes it impossible to use standard Wick's the- 
orem for constructing a diagram technique. Besides, the site-centered single-electron wave 
functions, used for the formation of many-electron ones, are not orthogonal to functions of 
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other sites. In the case of weak Coulomb interaction one does not require the use of a non- 
orthogonal basis. However, for many band structure methods this may still be desirable, 
for technical reasons. In this work we have shown that Schwinger's method^ of sources, 
which has become a standard method in solid state physics due to the work by Baym and 
KadanofiH can be used for constructing a regular perturbation theory for systems with 
strongly correlated electrons and we have shown that this leads to important advantages. 
On one hand, this approach allows one to overcome the difficulties existing in the original 
formulation of the diagram technique from the atomic limit@i! and to generalize the tech- 
nique for d- and /-electron systems. On the other hand, the approach happens to be fruitful 
also in overcoming the difficulties which come from the necessity to use a non-orthogonal 
basis set. We hope also that, it gives a convenient tool to derive effective low-energy equa- 
tions and to perform an analysis of corrections, arising beyond the standard Kohn-Sham 
equations. It is worth to note that since the theory is formulated in real space, it can be 
applied not only to the case of a regular lattice, but also to the case of impurity problems, 
or non-regular lattices. Of course, for the latter case an additional configurational averaging 
should be made. The perturbation theory developed here has a very close similarity to the 
ones for the boson and fermion systems. In particularly we have demonstrated this by the 
derivation of equations for the GFs: taking vertex in zero, first, and RPA approximations 
correspondingly, and we obtained the equations in the Hubbard-I, mean field, and RPA 
approximations. 

Some important questions, however, remain to be solved. It is a well-known fact, that a 
direct calculations of Hubbard's parameters for the Coulomb repulsion, for instance via the 
Slater integrals, give too large values. Within our theory the intra-atomic parameters enter 
the zero Hamiltonian, while the mixing and hopping become perturbations. The energies of 
intra atomic transitions of fermion type, which are determined by the bare value of Hubbard 
U in the atomic limit and in the Hubbard-I approximation, are renormalized by mixing and 
hopping in our theory already in the next, mean-field approximation. These renormalized 
parameters still determine the positions of the centers of the correlated bands and the 
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structure of the equations for the renormalized energies of transitions is almost coinciding 
with the one within well-known non-crossing approximation. It is interesting to note that 
mathematically the mechanism is analogous to the Lamb shift in electrodynamics, but here 
it is provided by conduction electrons. The energy difference between the centers of these 
bands can be interpreted as renormalized Hubbard U values. In a parallel work@ we have 
shown that, on one hand, the corrections, generated by the mixing interaction, can decrease 
or increase the value of effective Hubbard U depending on the particular solution chosen. 
On other hand, this mechanism gives small corrections and cannot provide the expected 
reduction of Hubbard U. 

In a series of works@El'Ei a separation of /-electrons into two subsystems, localized and 
delocalized, have been considered within the DFT approach on a phenomenological grounds. 
Application of the technique, developed in present work to this issueSi^l gives solid grounds 
to this assumption and provides also a better insight. Particularly, it is shown that the 
role of the localized /-wave functions is played by the combination (25) of the transitions 
[r n _x,r n ] with an effective single-electron wave function, 

^(iW„) (r) _ ¥ o Wer{r) s J d(p 1 )d( ra )...d( ril _ 1 )«rj w _ i ((p lj ra , r^^r.CCn, r a , .... r n _ a , r), 

(123) 

corresponding to the lower-energy transition A 1; while the delocalized /-wave functions 
correspond to combinations of wave functions 
y(r n ,r n+1 ) {r) = ^u PP er {r) = I d ( ri ) d ( ra )... d ( rn )^ n (( ri) r2) r n )* r „ +1 (( ri , r 2 , r n , r), 

(124) 

The question how much the intra-atomic Coulomb / — /-interactions are renormalized by 
those corrections from the (spd — /)-interactions, which are not included in the standard 
LDA scheme is still open. One can expect that these interactions should provide a strong 
reduction of the effective U*-s. This question is especially important for a calculation of 
the self-interaction correction (SIC), which plays a decisive role in the problem of volume- 
collapse transitions in rare earth metals^ and the theory of Mott insulators. We expect that 
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the expression for this correction, which arises in the equations of motion for the GFs when 
the self-energy is approximated by the difference between the static part of RPA-screened 
Coulomb interaction and the matrix element of the LDA potential, should be accurate. No 
problem with double counting and orthogonalization arises within such an approach. 

The method suggested here is a special type of perturbation theory, which does not 
contain any approximations in its general formulation. The Hubbard-I approximation (HIA) 
and mean field approximation from atomic limit (AL-MFA) can be made with reasonably 
moderate efforts. It is relevant to emphasize here the similarities and differences between 
our AL-MFA from the popular method of dynamic mean field (DMF)i. The AL-MFA 
differs from the HIA by taking into account the renormalization of energies of intra-atomic 
transitions by mixing interaction and hopping. Therefore, besides self-consistent equations 
for many-electron population numbers, which arise in HIA, the system of equations for self- 
consistent centers of gravity of the Hubbard sub bands is present in the calculation in AL- 
mfaH. A renormalization of this type is accounted for also in the DMF theory. However, it 
is not clear at the moment, to what extent they are different. Eq.(122) for the renormalized 
energies of intra atomic Fermi-type transitions in the AL-MFA is non-linear and sometimes 
contains more than one solution; the number of poles in the fermionic GF (or, number 
of Hubbard sub bands) is equal to the number of intra atomic fermionic transitions only 
near atomic limit. The DMF method is analogous to the coherent potential approximation 
and, therefore, can give some other poles. To be specific, the DMF gives a central peak 
in the density of states of fermionic excitations for the paramagnetic state of the s-band 
Hubbard model. The AL-MFA equation for renormalized energy of transition also may 
have second solution in vicinity of Fermi energy!!, what can provide the central peak in the 
density of electron states. However, in many cases the upper transition A 2 being slightly 
above Fermi energy hybridizes with non-/ bands and gives birth to some tails under Fermi 
energy. These peaks are seen in electron spectroscopy^ and are close enough to Fermi energy 
in order to mask weak Kondo peak. The details of the mechanism of formation of these 
additional bands in the region of moderate coupling constants in terms of Hubbard-band 
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excitations remain to be elucidated. There is a difference between DMF and our approach in 
the treatment of spatial correlations. Namely, in applications of the DMF approach only a 
local, momentum- independent, self-energy has been considered, whereas in our theory first 
(k, cu)-dependent corrections in the expansion of the self-energy arises in the RPA (due to 
kinematic interactions). 

At last, we note that it is easily seen from the close similarity between our construction 
of perturbation theory and the one by Baym-KadanoflS that the description of interacting 
delocalized and localized electrons in non-equilibrium states needs only slight generalization 
of our theory. Example of calculation of effects of non-orthogonality on the tunneling pro- 
cesses via a region with strongly interacting electrons is given in Ref. 48. An application of 
the technique to spin systems will be given elsewhere. 
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APPENDIX A: PARAMETERS OF THE HUBBARD- ANDERSON MODEL 

Let us now discuss the parameters of the model. The energies E® according to defini- 
tion, are E® = (T\H ion \T). For the case of Russel-Saunders coupling, i.e. for rare earths, 
corresponding expressions are given by Eqs. (7. 18), (7. 19) in RefM Main contributions to 
these energies come from pure Coulomb interaction, 

5E V , C = \n{n - l)F® - (Y\v LDA \Y). (Al) 
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Typically, the bare are very large, ~ 20eV, what lead HubbardB to the conclusion 
that only the states with n = n ,n Q ± 1 from the parabolic dependence of -Er( n ) can give 
contribution to the Hamiltonian. Here no is dictated by the valence of the ion. As has 
been explained by McMahan0, the interaction with non-/ electrons strongly decreases the 
Umodei = (Ar n+1 ,r„ - (Ar n ,r„_i), where A r>ri = (E r - E Tl ). It worth to emphasize the 
difference between our and Hubbard's schemes. Hubbard and, also, Irkhin, diagonalized all 
on-site electrons, treating all of them in equal footing; this means that they have accepted 
the assumption that all on-site electrons are in strong-coupling regime. We assume that 
the conduction electrons (which are above tops of Coulomb barriers, can be described as a 
weakly-coupled subsystem. 

Let us now discuss the mixing interaction. The matrix element W in Eq.20 consists of 
three terms: W mix = W {1) + W {2) + W {3) . The term W {1) ~ Hc + f + H.c. describes the 
single-electron LDA-type of mixing, which is present in Eq.(8). This matrix element is: 

W jL,j'a> = S ^HjLji m > l a{f m > l a) a = ^ H jL j'm\o (^n \ fm\a I ^n+l) ( A2 ) 

The other contribution to B. m ix arises from those terms of (H — Hlda) which have the fol- 
lowing operator structure; c + X a and X a c. Obviously, the terms of the Coulomb interaction 
8H%* = C{U: ft ft he, + ClMftftc.h + C2U 1 f ciftf 2 fi + CMlftcthh, which con- 
tain the /—operators related to the same site, should not be decoupled. We rewrite these 
products of /-operators in terms of X-operators, 

Utfth) = HUtfthfx\ (/ 3 + / 2 /i) = E(/ 3 + / 2 /i) a x a . (A3) 

a a 

In this place a possibility exists to consider a simplified models, omitting some terms of 
genius Coulomb interaction, but keeping corresponding terms in LDA. The Hartree part of 
the LDA potential contains these terms in the decoupled form. If we perform the decoupling 
of these terms in the Hartree- Fock fashion, there appears also exchange terms. However, the 
exchange interaction is treated within LSDA theory by the term ~ (ap 1 ^ 3 (r))j L j L iat L ajLi in 
the potential. Therefore, it is reasonable to regroup the terms in order to represent every 
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term of the perturbation Hamiltonian in the form of a deviation of the non-approximated 
terms from the ones used in the LDA approach, as has been discussed in Sec. II. In order to 
avoid double counting we should not touch other terms of the Coulomb interaction if we do 
not consider its LDA-counterpart, if we do not take into account some of the Coulomb terms 
(if — H LDA )^° U \ corresponding terms enter all equations in the LDA form. Particularly, 
we treat in this way the terms of the form ~ Cc + c + cf and similar ones. Then we come to 
the following form for the remaining terms of the mixing matrix elements. For the density- 
density part we have, 

W (2) = \ - r cfff 

jLJ'a' ^jL,j'm ll o J ;j'm l 20 J ;j'm l3 <r 
{m ; },<r 



X Kfmna'fmuv'frnisa) {fmi 3 a') (fmna'fmnv' 



LDA] 



~ V c cfff 

Z-^i jL,j'maa';j'mi2(T';j'mi3a 
{m ; },<7 

], (A4) 



and for the exchange part, 



(3) _ \ - r cfff r X Jf W LDA 

iL,j'a' ^jLj'miia'-j'mi2a';j'm l3 a u m2,rni3 u <T<T\Jm n a' ) "'j'mna' 

{mi}/ 



-Y,Mr)) jL , mi Af mil „T'- (AS) 

As seen the term ~ W^ j/a , = for the empty orbitals. It is this matrix element that is 
responsible for the different mixing of conduction electrons with the intra atomic transitions 
/ n_1 — > f n and f n — > f n+1 . Obviously, this is a pure correlation effect: decoupling within 
the weak-coupling perturbation theory will make it vanish. 

The prohibition to decouple the single-site products of /-operators leads to Coulomb 
renormalization of the LDA-hopping too: 

tjj' =t tii' + t 2jj'i ( A6 ) 

tfjf = Efe' - (vLDA)pM) a (fu) b , (A7) 
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t: 



ab 
2jj 




[V 



(A8) 



Note the t2-contribution from Coulomb interaction works similar to corresponding contribu- 
tions to mixing: it makes hopping very different for, say, transitions [r„_i, r„] and [T n , T n+ i]. 

APPENDIX B: EXAMPLE: THREE-ELECTRON STATES AND 
COMMUTATION RELATIONS 

Let us consider a simple example of three electron states and corresponding commutation 
relations for the Hubbard operators, constructed on these states and belonging to the same 
site. A general three-electron state |T) can be written as follows: 



Here ^-functions are needed in order to establish a sign convention, i.e. in which order have 
to be taken the /-operators in the many-electron wave function; 9 V1V2 — 1 if v\ > v 2 and 
zero otherwise. Then, according to the definition of the X-operators, X rr ' = |r)(r'|, the 
multiplication rule is: 




rTi -^rar' 



(r\|r 2 )x 



(B2) 



From (B.l) we have 




(B3) 



Then, from the anticommutation relation 




i-i 

71)1,1711/ 1 



(B4) 



For the same site n = m and we see that 



= det 




(B5) 
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and, therefore, the commutation relation is given by 

{x rr 1)X r 2 r' } = ^ + ^x^. (B6) 

The equation of motion for the operator Xf Fl contains the commutator 

£ r[Xf lF2 , h^ r ] = J2[vir 2 ,nr£rX[n r - ^r^r^rA^ 2 ], (B7) 
n,r n,r 

where 

Af/^r^ri. (B8) 

According to our assumption the direct overlap integral on wave functions of /-electrons 
belonging to different sites is negligible, therefore, main contribution comes from the terms 
with i = n, since i^r 2 ,nr involves determinant of the inverse of overlap matrix on / — /- 
matrix elements and each of them is proportional at least to second order of overlap between 
wave functions of band electrons and /-electrons. Without this approximation one has to 
solve problem of interacting via conduction electrons clusters. In this work we use the 
approximation 

J2 er[Xf ir *,h™] « Y>^evX^ v - w^X?*] » A r2ri X rir2 (B9) 
n,r r 

with 

Ar 2 r\ = ^ir 2 ,iT2 £ r 2 ~ ^ri.ir^ri. (BIO) 



APPENDIX C: EQUATIONS OF MOTION FOR OPERATORS 

1. Equations of motion for quasi-fermion operators 

In order to construct equations for the Green's functions we need the equations of motion 
id t F = [F, H — fxN] for (quasi)fermions. As seen from Eq.(28), the particle number operator 
has the same operator structure as the Hamiltonian (26-27). Therefore, in order to obtain 
H = H — fiN, we should make the simply the following shifts: 
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E° - E° = E° r -fi^MMh, (CI) 

M 



HjL,j'L> -> H jL j /v = H jL ji V - fiO jLJ > L >, (C2) 

""^If =t ii' _ ^ E ( 1 _ S jj')(fM 1 ) a OjM 1 ,j'M 2 (fM 2 ) a ) ( C3 ) 
MiM 2 

and 

WjLJ'a -> WWa = W^-LJ'a ~ A* E <Wm(/m)", (C4) 

M 

in the Hamiltonian if. Below we assume that such a shift has been made, and will accord- 
ingly drop the tilde. Then the equations of motion for the c- and X-operators can be written 
as: 

idtCjL = E (0- l H) 3Unh2 c nL2 + (M c ) jL , (C5) 

i^x; = a s x; + (m x ); (ce) 

where for a = [r n _i, T n ] the energy of the transition A s = -Er„^r n ,r„ - ^rw^^rw- The 
coefficient v — 1 when the /-orbitals are orthogonal to other and v < 1 when they are 
not. The operator M describes deviations from the solution at zero mixing interaction in 
the particular potential chosen (say, LDA). The M-operator for the c-subsystem consists of 
linear and non-linear parts, 

(M c ) jL = (M l c in )jL + [{Ml) jL - (Mj) jX ], (C7) 

where the linear one is given by 

(M i : n )jL=j:w jL , j2L2 x%, (cs) 
J2L2 

and 

WiLj*» = {0~ l W) jL , j2a2 + E (0- l ) jL , 2M {f M ) a ef 2 El T . (C9) 

r,M,a 
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The non-linear terms are given by the operators (M})j L , (M c )j L which can be written in the 
form 

(m}) jL = Eio-'UjMfMTP^t?^ + w; iaj2L2 c nL2 ] (cio) 

and 

(Mj) jX = DO-W(/ M ) c [W + W nLl , na c] lLl ]P™. (Cll) 

Here P a ' S = {X a ,X 5 } = ea,6 ? Z«, P a ' b = {X a ,X 6 } = €a\Z*. The expression for the 
M-operator of the /-subsystem can also be written in the form 

(Mxh = [(Mi); - (M±)J] (C12) 

with 

(^4)1 = E + E ">',,,/,- 7 7''':/:-/,, (ci3) 
(m x T 3 = E*£4^ CQ • E^/,,,^/,/7' 

Here the matrix elements are given by the following combinations 

< b , j2L2 = w; ibJ2L2 + ^Uti)\o-^H) mnL2 , (ci4) 

M 

tt = 4 + UfD\0- u W) jM , 2C (CIS) 

M 

u,i = («,T)t > f i = ( f T)t. (C16) 

Some features of the correlated electron system can be seen already from these equations. 

First, we see that instead of equations for /-operators, as in case of non-correlated 
electrons, we have to write equations for their complicate combinations X a . This leads 
to an unpleasant consequence: the number of variables and, correspondently, equations, 
increases drastically: every /-operator is split into many Hubbard sub-amplitudes. The 
number of transitions involved depends on the valence of the /-ion. 

Second , we find that due to non-orthogonality the /-level also contributes to the mixing 
interaction: the term in Eq.45 describing linear mixing would disappear for the orthogonal 
basis. 
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Third, the non-orthogonality transfers non-linearity into the conduction electron subsys- 
tem. This correctly reflects the underlying physics; the c-electron spends part of its life as an 
/-electron, therefore, the mixing and hopping which it experiences should be correlated to 
the same extent. It is this reasoning which forced us to perform operations in this particular 
order. First, we take into account correlations, and then perform orthogonalization. The 
other motivation is of a technical nature; this construction allows us to avoid consideration 
of a combinatorics of tails for vertexes with one c-tail and many /-tails. On the other hand 
if we now decide to write the equations in terms of /-electrons (i.e., fermion /-operators), 
some features of the many-tail vertexes are clarified by these equations. Particularly, the 
many-electron GFs, which are not allowed to be decoupled, can be obtained simply by 
transformation of all Hubbard's operators back to the /-operator form. 

Fourth, the equation for the X-operator does not contain linear mixing. 

Fifth. The right-hand side of the equation has some similarity to the Boltzmann equa- 
tion; both contributions, from mixing and from hopping, consist of two terms, "incoming" 
and "outgoing"; we denote them by the superscripts j and j, correspondingly, in order to 
emphasize that the first describe intra-atomic transitions " up" , from the states with n elec- 
trons to the states with n + 1 electrons; the j-terms describe the inverse processes. As we 
have seen above (see, e.g., discussion after Eq.110) later, the contributions from the latter 
terms always contain a small parameter. 

2. Equations of motion for quasi-boson operators 

Let us introduce the vector-operator rj, which will describe all fermion or fermion-like 
transitions in our system. For example, for the system with n conduction bands and m 
intra-atomic fermi-like transitions it is given by: 

Vj\ = ( c ]lv •••' c iLn'^i 1 ' •••>^j m > c j£i' •••'S'in'^i 1 ' x j m )- ( C17 ) 
As is seen from this notation, for r] = c or the index A takes a value L while for r\ = X 
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or X\ and A can either be b or The Hamiltonian without external fields generates the 
following equation of motion for the operators ; 

id t Z] = £4rA^} + Vi Lujb cl Ll X» + Vi niL X]c n2L2 
r 

+^1101,71202 [^"iJ^ai,6 + ^n2,j^a2,b]^nj^n2' 
= E + K?iAi,n 2 A2^iAi^2A2- (C18) 

r 

Let us remind the reader that for the diagonal operators (i.e., for the transition £ = [r, T]) 
the first term in the right-hand side of Eqn.C18 vanishes, since for this transition A^ = [ r r j = 
u rr E r — u rr E r = 0. Furthermore, we have defined the constants of the interaction in such a 
way that V^ lXl naA2 = when (niAi) = {n\L\) and (n 2 A 2 ) = (n 2 L 2 ) simultaneously, i.e. when 
both indices relate to c-electrons. In all other cases they denote the following combinations 
of the matrix elements : 



A r = 

]b,n 2 L 2 


OjM,n 1 L 1 U'M) C£ - b ' i h niLun2L2 + W* aun2L2 el ai , 


(C19) 


V ( T b = 


~^n 2 \ 2 ,jM^ MT £ b'^ l nxLx,n 2 L 2 + ^V ni L 1 ,ja 2 £ b 2 > 


(C20) 


V*- h = 


— O^ 1 (f Yr c ^W* 

n 2 L 2 ,jM\J M / °b vy n\a\,n 2 L 2 i 


(C21) 


v<- h = 


~^jL 1 ,n 2 a 2 C ) jM,n 1 L 1 (fM) C£ l^- 


(C22) 



APPENDIX D: AUXILIARY FIELDS A JXJ , X ,(T) 

The auxiliary fields Ajx,j'y(t) have the following non-zero elements: 
For r/x = c we have 

A< fSa(t) = E £b a%^)UM)\0-% L , 3lM (Dl) 

j'Mb,^= n ,A,r 

and 

A< fL,Pa(t) = E £b ^fMfM) b (0- l ) jL , fM . (D2) 
j'Mb,fj 

For r/x = X we have 
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4S3(*)= E ^{1)5^ = 5^% (D3) 

and 

= E^m? = ^- b . (D4) 

For r^A = we have 

= - E 4^(0 (/mHo- 1 *)^ (ds) 

and 

= - E ^-(oa^^o-^w. (D6) 

ryMa 

For 77* = X" 1 " we have 
and 

45? W = S ir E (*) = 8 jr A%. (D8) 
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FIGURES 

FIG. 1. Corrections to the self-energy in the approximation Hubbard-I (Eq. 78). The dashed 
wiggle line denotes the part of mixing 6 r , the triangle with inscribed a denotes the a-matrix, given 
by Eq.64, the mixing interaction is denoted by the wiggle line (see Eq.84). The expectation value 
(Zs b ) is denoted by circle-cross on the top of the a-triangle. Summation over bosonic index If, is 
implied in this top what transforms the product ai2,3 b (Z) into P\i- 

FIG. 2. Simplification of the graphical notations in the case of equations for fermions. In the 
vertex containing a-matrix the line of mixing interaction goes out of the vertex, while in the one 
containing /3-matrix it enters the vertex; /3-matrix appears always together with the r-matrix (see 
Eqs.64, 72 and 88), which inverts direction of the line of mixing. This vertex describes the process 
involving two-electron (r n _2 — > r ra )-transition. 

FIG. 3. i?'-rule of taking functional derivatives of the pseudolocator line; changing a by (3 gives 
iT-rule (see Eqs.85 and 86,87). 

FIG. 4. Mean field corrections to the self-energy of pseudolocator. Note that both vertexes 
have three tails of pseudolocators and one of interaction; left vertex has the outgoing interac- 
tion tail, while right one has incoming interaction tail. The point from which the dashed line of 
r/-pseudolocator and the wavy line of the n — ^-interaction V12 go out denotes the r-matrix. 

FIG. 5. The second-order corrections to the vertexes from kinematic interactions. 

FIG. 6. Corrections to the self-energy of the ??-pseudolocator generated by the second-order 
corrections to vertex. 
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